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vector of gravitational forces 

angle between planes formed by unit vectors from spacecraft to 
planet and two stars 

bias error 

declination of planet 

expected value 

orbital eccentricity 

orbital angular momentum 

unit vectors along Cartesian axis 

maximum ratio of interpolation error to measurement error 

right ascension of planet defined with respect to star-fixed 
coordinates 

range vector from spacecraft to center of planet 

magnitude of r 

half the planet subtense angle 
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unit vector along r 
Cartesian components of u 
velocity vector 
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PRELIMINARY DETERMINATION OF SPACECRAFT TRAJECTORIES 


USING ON-BOARD OPTICAL MEASUREMENTS 
John D. McLean 
Ames Research Center 


SUMMARY 


This report deals with the problem o£ determining the trajectory of a 
spacecraft from on-board optical measurements without a priori information. 

Only sufficient information to determine the trajectory is used, and it is 
assumed that this preliminary estimate will be improved later by the use of 
differential corrections and redundant data. 

Five basic types of measurements, representing instruments considered 
feasible or within the state of the art, are considered. The effects of ran- 
dom and bias type measurement errors and interpolation errors arising from 
nonsimultaneous measurements are analyzed. Only observations of the central 
body and two stars are considered. 

Within certain limitations all five measurement types will give 
satisfactory results. Two important factors are the proper choice of stars 
and the calibration for eliminating bias errors. The errors for four types of 
measurements are strongly dependent on the relative location of the projection 
of the radius vector from the vehicle to the central body onto the plane 
formed by the unit vectors directed from the central body toward the two stars. 
The choice of stars is greatly enhanced by allowing measured angles as large 
as 90°, substantially larger than the present state of the art. Irradiance is 
a major source of bias error, and every effort should be made to correct for 
it; even a crude correction would be helpful. 

The random position error can be approximated by the radial error from a 
single subtense angle measurement, and the accuracy of this approximation 
increases with range from the central body. This approximation can be used to 
show that within the present state of the art, the maximum distance from the 
central body at which the measurements considered in the study would be useful 
is about 200 times its radius. The usefulness of the methods at ranges of 
less than /I times the radius of the central body, while theoretically 
possible, is of questionable practicality. 


INTRODUCTION 


Many navigation schemes require an initial estimate of the trajectory 
and use very sophisticated data processing techniques. Such systems are sub- 
ject to partial failures, which, if not discovered immediately, make the 
information in the on-board computer useless. It is thus essential to have 
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available a simple back-up navigation scheme that can determine an orbit to 
a reasonable accuracy with no a priori information. 

So far, studies of back-up navigation methods have concentrated on the 
use of manual computations and on the calculation of perigee altitude (e.g., 
refs. 1 and 2). A more desirable system, however, would consider all elements 
of the trajectory. Reference 3 describes such a system, which is designed for 
near-earth orbits of very small eccentricity. This approach uses manual cal- 
culations and a large amount of precomputed graphical data. For more general 
applications, such as lunar and interplanetary trajectories, greater accuracy, 
in terms of significant figures, is required, and it appears impractical to 
obtain this accuracy without the use of redundant data. 

The approach used in this study is the same as that used by astronomers 
in determining the orbit of a newly discovered body. That is, a preliminary 
estimate of the trajectory is made using only enough data to determine the 
position and velocity vectors at a given time. The estimate is then improved 
by the use of redundant data from additional observations. This technique 
does not require any greater accuracy in measurements but it does require the 
use of a small automatic computer. 

This report is concerned with the first part of the problem, that is, 
with a preliminary determination of the trajectory, or, in other words, the 
position and velocity vectors of the spacecraft at some given time. The 
determination of velocity consists essentially in making two position fixes 
and then interpolating to find position and velocity at some intermediate 
time. 


Hamer and Mayo (ref. 4) reported on the accuracy of obtaining position 
fixes using a number of different combinations of simultaneous on-board mea- 
surements including radar ranging. This report presents a more detailed error 
analysis of several possible methods that require fewer measurements and 
include some practical constraints. This analysis is carried out by the 
linear perturbation method, which Hamer and Mayo have shown to be valid for 
fairly large errors in the type of measurements considered. 

Only optical observations of two stars and the central attracting body 
involving measured angles no greater than 90° are considered. Five different 
types of measurements are investigated. The instruments used are either 
within the present state of the art or are considered to be potentially feasi- 
ble. In some cases a complete position fix cannot be made with a single 
observation, and it is then necessary to refer the data to a common epoch by 
interpolation. The report considers errors due to interpolation as well as to 
measurement inaccuracy. 

Since the velocity is determined directly from position-fix data, most of 
the error analysis deals with the accuracy of position fixes, and the errors 
in determining the corresponding velocities require only a brief treatment. 

The object of the error analysis is to determine whether the preliminary esti- 
mate of the trajectory is accurate enough to allow the use of linear perturba- 
tion theory to process subsequent redundant data for improving the estimate. 
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It is difficult to assign numbers to such an accuracy criterion since the 
errors in the linear perturbation theory depend upon the time before addi- 
tional data are acquired and the nature of the trajectory as well as the 
errors in the preliminary estimate. In most cases, however, the preliminary 
estimate should be accurate enough if the magnitude of the errors in position 
and velocity do not exceed 2 or 3 percent of the magnitudes of the total 
position and velocity, respectively. 


DESCRIPTION OF MEASUREMENTS 


The measuring instruments for a back-up navigation system should be 
simple and reliable; therefore, only optical instruments have been considered. 
In addition, only observations of two stars and the central body were used in 
order to simplify the measurements and computations. Triangulation of two 
celestial bodies, which usually gives the range more accurately than the cen- 
tral body subtense angle (ref. 4), was excluded. This was done because of the 
number of measurements required and the difficulty of interpolation if the 
necessary angles cannot be measured simultaneously. 

The error analysis methods described later in the report are valid for 
angles as large as 180° between the two stars and between each star and the 
center of the central body. For the computation of data, however, these 
angles were restricted to magnitudes of 90° or less, giving 90° total field of 
view. This restriction is reasonable since it can be shown that if stars are 
available anywhere within this allowable measurement range, then accuracy will 
not be improved by increasing this measurement range. Angles up to 90° are 
included for completeness, even though this limit is beyond the present state 
of the art. For example, the Apollo and Gemini spacecraft windows restrict 
the angles measured with the hand-held sextant (ref. 5) , to be within a field 
of view of about 50°. The Apollo space sextant can measure angles up to 57° 
from the principal axis to give a field of measurement of 114°. Even so, the 
fact that the measured angles cannot exceed 57° prevents the use of many pairs 
of stars within the field of view. 


Definition of Variables 

We wish to determine the position of the spacecraft in a Cartesian 
coordinate system having the origin at the center of the dominant attracting 
body. However, for convenience we will translate the origin from the central 
body to the spacecraft as shown in figure 1. The vector f in the figure is 
the position vector of the planet^ with respect to the spacecraft and is the 
negative of the desired position vector. 


^In the remainder of the report the word "planet" will be used to refer 
to the central attracting body even though that body might actually be the 
Sun or Moon. 
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For convenience the coordinate 
system used will be a Cartesian sys- 
tem as follows: The positive X 

axis lies in the direction of star 
number 1, while the Y axis is in 
the plane of the two stars, normal 
to the X axis and positive when on 
the same side of the X axis as 
star number 2. The Z axis com- 
pletes a right-handed orthogonal 
system. This report deals only with 
finding r in the coordinate system, 
but the coordinates can be trans- 
formed readily into any inertial 
system in which the directions of 
the stars are known. 


The angles RA and D in 
figure 1 are, respectively, the 

Figure 1.- Geometry tor position fix. right ascension of the planet mea- 

sured in the X-Y plane from the 

X axis and its declination above (or below) the X-Y plane. It is pointed 
out in reference 4 that another star must be observed to avoid ambiguity. The 
ambiguity occurs here as an uncertainty in the sign of D, but the observer 
can determine the proper sign easily from the relative positions of the planet 
and the two stars. (If one draws an arrow from star 1 to star 2, D is 
positive when the planet is to the left of the arrow.) 


The "star- center" angles ijii and i(i 2 are measured from stars 1 and 2, 
respectively, to the center of the planet, while S is defined as half the 
planet subtense angle. 


We define the corresponding star-limb angles Yi and respectively, so 

that 


ijii = Yi + kiS 


( 1 ) 


where 


ki = ±1 

The positive sign is used for Yi measured to the limb nearest the star and 
the negative sign for the far limb. The angle between the planes formed by 
r and the lines of sight to the two stars is denoted by B, and is the 
angle between the two stars. 

If u is a unit vector in the direction of r then 


r = ru = r 



( 2 ) 
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and 


uj = cos RA cos D 
U 2 = sin RA cos D 
U 3 = sin D 


( 3 ) 


while 


r = R CSC S 

where R is the radius of the central body. If we normalize by using the 
planet radius as the unit of length, then R = 1 and 

r = CSC S (4) 


Measuring Instruments 

Three general types of measuring instruments were considered. The first 
type would measure RA, D, and S or ijji , 3 simultaneously. The 

photographic technique investigated by Walsh (ref. 6) could be used for this 
purpose. Walsh's approach uses a sextant-type measuring instrument with the 
human eye replaced by photographic film. Sketch (a) illustrates one possible 
example of such an instrument for two dimensions. The star and central body 
need only be within a relatively large field of view, and the incremental 
angle, which must be added to a (the angle between the optical axes), can be 
measured from the film. The subtense angle and center of the planet can be 
determined from the film and, therefore, the star-center angle can be measured 
directly. Walsh's results indicate the photographic method to be potentially 
feasible, and although he considered no more than two lines of sight, the 
method could be extended to three. The three- line-of-sight method would allow 
the simultaneous measurement of S and either and 4^2 3-nd D, 

depending on the configuration of the instrument. The latter combination. 
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which could also be measured by a theodolite mounted on a stable platform, 
will be used as a standard fof comparison in this study. 


Field of view of principal 
line of Sight 


The second type of instrument considered is capable of measuring two or 
three angles simultaneously without the use of photography. The trisextant 
reported on by Novak (ref. 7) is an example of such an instrument. The tri- 
sextant has three lines of sight and the principal line of sight has (or could 
have) a field of view large enough to observe a substantial portion of the 
planet, while the two secondary lines of sight observe stars. The operator 
sets the images of the stars tangent to the disk of the planet simultaneously 

as shown in sketch (b) . The readout 
scales on the instrument provide the 
angles Yi> Y 2 j B. The portion 

of the planet between the two stars 
must be within the field of view of 
the principal line of sight, and this 
requirement restricts the possible 
location of stars when the spacecraft 
is near the planet. This restriction 
may be dealt with in a number of 
different ways, such as using a 
larger field of view and less magni- 
fication at shorter ranges, and has 
not been considered in this study. 

An experimental model of the Martin 
trisextant has been tested at Ames 
Research Center and found to be fea- 
sible from the operational standpoint. However, this model was not suffici- 
ently accurate for space navigation, and the ultimate accuracy obtainable with 
such an instrument is unknown. 



The third instrument considered, which measures only a single angle at a 
time, is typified by the hand-held sextant described in reference 5. Data 
obtained from the instrument in ground tests and space flight (refs. 5 and 8) 
were found to have random errors of less than 10 arc seconds standard devia- 
tion. Bias errors ranged in magnitude from 0 to 28.7 arc seconds, depending 
on the measurement conditions. Further studies may lead to the reduction or 
calibration of part of the bias errors. 

In tests of instruments such as the hand-held sextant it has been found 
difficult to estimate the planet center accurately, especially when the full 
disk is not visible. Therefore star-center measurements are considered only 
when photography is used. 


Types of Measurements 

In order to determine the position vector, at least three of the angles 
defined must be known at a given time. The consideration of different practi- 
cal combinations led to the following five types of measurements. 
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Type 1- Simultaneous measurement of the right ascension, RA, 
declination, D, and subtense angle, 2S, by means of photography (or an 
inertially stabilized theodolite}. As will be seen in the section on Error 
Analysis Methods, this is the only one of the five types of measurements for 
which the error equations do not contain singularities. Therefore type 1 
measurements have been used as a standard of comparison for the other types. 

Type 2- Simultaneous measurement of yi> Y2» 2s. This measurement 

would require photography, but it could be approximated by using the trisex- 
tant to measure 2S, then Yi and Y2 simultaneously, followed by another mea- 
surement of 2S. Linear interpolation, discussed later in the report, would 
be used to refer S to the same epoch as yi and y 2- This approximation is 
not analyzed. 

Type 3- Simultaneous measurements of yij Y2> and angle B, between the 
planes of the two star- limb angles. This single measurement with the 
trisextant would allow a complete position fix. 

Type 4- Simultaneous measurement of one star- limb angle and the subtense 
angle followed by simultaneous measurement of the other star-limb angle and 
the subtense angle at a later time. Interpolation is used to refer these two 
sets of measurements to a common epoch. This is another possible application 
of the trisextant, or it could be done with the two-lines-of-sight 
photographic sextant described by Walsh. 

Type 6- Measurements of the two star- limb angles and the subtense angle, 
one at a time in a suitable sequence, followed by interpolation to refer the 
measurements to a common epoch. These measurements, which could be performed 
with the hand-held space sextant described in reference 5, constitute the only 
one of the five types that is truly within the present state of the art. 


DETERMINATION OF POSITION AND VELOCITY 


This section of the report presents the equations used for determining 
position and velocity from the measured angles. The angles RA and D are 
measured directly only in type 1 measurements, and for the other types the 
ui in equation (3} must be expressed in terms of the angles (yi, Y2> ^ or B) 
actually measured. However, it is more convenient to combine the Yi with S 
using equation (1) in order to get the star-center angles The resulting 

equations for the uj^ and the equations for determining S from type 3 mea- 
surements are derived in appendix A. Types 4 and 5 also require interpolation 
formulas for determining position and these are derived in appendix B along 
with those used for determining velocity. 

The position-fix equations for all types of measurements are summarized 
in the following paragraphs after which the computation of velocity is 
discussed. 
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Position Fix Equations 

Type 1 measio’ements- In this case equation (3) can be used for u while 
r is given by equation (4) . 

Type 2 measiai’ements- Again, r is given by equation (4), but u is given 
by equation (A8) as 


cos ij;^ 

cos \p 2 - cos cos (p 

^ sin <p ^ 

± Jl - U^^ - U2^ 

The third component is equal to sin D and the sign is chosen by the method 
described earlier. 

Type 3 measurements- For these observations S is computed in terms of 
Yij Y 2 > B using the following equation from appendix A; 


cos(2S + k^Yi + 


2 cos (j) - cosfkjYi ■ k 2 Y 2 ) (1 ^ik 2 cos B) 
1 - kik 2 cos B 


( 6 ) 


where k^ and k 2 are defined by equation (1). Once S is known, and ip2 
can be found using equation (1), and u is determined by equation (5). Equa- 
tion (A15) is an alternate equation for u 


cos 

'Pi 


cos 

'Pz 

- cos iJjj cos (j) 



sin (j) 

±sin 

'Pi 

sin Tp 2 sin B 



sin (j) 


Type 4 measurements- This is the first case in which it is necessary to 
use interpolation. Measurements of one star-limb angle and the subtense angle 
at each of the times are indicated in table 1. 


TABLE 1.- REGULAR SCHEDULE FOR TYPE 4 MEASUREMENTS 


Time 


^2 

^3 

Angles measuijed 

Y ^ and S 

Y 2 and S 

Y and S 
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Equation (4) can be used to compute r at each time in the table while 
equation (1) gives at ti and t 3 and ip 2 at t 2 - 

Next Gibb's method (ref. 9) is used for interpolation to find (t 2 ) . 

The interpolation formula, derived in appendix B, is equivalent to third order 
in time. For convenience t 2 is assumed to be zero, in which case 

1 r - - A(t,) 

r(t2) = tarCtJ - tir(t3) ^ — titsCti - 2t3) 

A(t3) 1 

^ tit3(2ti - t3) (7) 

where A(t) is the gravitational force acting on the spacecraft at time t. 

At times ti and t 3 , r and have been determined so that if w^ is 
the unit vector in the direction of star number 1 it is possible to compute 


r(ti) • wi = r(ti)cos (ti) 


Thus equation (7) can be used to find cos 4^i(t2) by taking the dot product of 
both sides with w^ if a conic trajectory is assumed. In this case 


A(t) 


pr(t) 

r3(t) 


and, using equation (7), 


A(t) • wi = r-— cos ij;i (t) 

r^(t) 


This approximation is necessary since all three components of r(t) would be 
needed to compute A(t) • w^ for the n-body case. The resulting expression 
for cos ipi (t 2 ) is 


cos >J^l(t 2 ) 


1 

ts - tl 


r(t3) 

r(t2) 


(r(ti) r Ftit3(ti - 
lr(t2) [3 ■ 6r3(ti) 

r ^ yt^tg (tg - 2t^) 1 

L ^ 6r3(tg) J 


2t3) 


cos 


■J cos ij;i (tl) 


( 8 ) 


Since cos ip 2 (t 2 ) and r(t 2 ) can be computed directly from measurements at t 2 , 
enough information is available to compute r(t 2 ) using equations (2), (4), 
and (5) . 
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Type 5 measurements- In this case it is possible to measure only one 
angle at a time, and it is necessary to depend partly on linear interpolation. 
Either the range, r, or one or both of the star- limb angles may be interpo- 
lated, but it will be shown later that linear interpolation is somewhat more 
accurate for r than for y- Therefore, from the standpoint of interpolation 
accuracy, the schedule in table 2 was considered best. 


TABLE 2.- REGULAR SCHEDULE FOR TYPE 5 MEASUREMENTS 


Time 


^2 

^3 

t4 

ts 



Angle measured 

S 

yi 

S 

^2 

S 


S 


Linear interpolation gives the range at t 2 , tt^, and tg as 


r(ti) 


(tj+i - ti)r(tj-i) + (tj - ti-i)r(tj+i) 

(^i+l) ■ (ti-i) 


C9) 


Once r(t 2 ), r(t 4 ), and r(tg) are known, equation (8), with the proper sub- 
stitutions for the ti, can be used to compute cos (t 4 ) . Then cos 
is determined from a direct measurement, and equations (2), (4), and (5) give 
r(t4) . 


The procedure just outlined uses linear interpolation to determine r at 
times t 2 , t 4 , and tg, after which Gibb's method is used for interpolation to 
determine cos if)i(t 4 ). Two alternate schedules which use only linear inter- 
polation are shown in table 3. The same time indices are used for comparison 
with table 2. For case 1 equation (9) is used to compute r(t 4 ); yCt) is 
substituted for r(t) in equation (9) to obtain Y;^(ttj) which can be used with 
r(t 4 ) and Y 2 (t 4 ) in equation (1) to give (t 4 ) and case 2 equa- 

tion (9) is used to compute Yi (^ 4 ) and Y 2 (^ 4 ) » and S[t 4 ) is known from 
direct measurement. Each of these alternate schedules uses equation (9) twice, 
whereas the regular type 5 schedule uses equation (9) three times and equa- 
tion (8) (from Gibb's method) once. Therefore, the alternate schedules 
require less computation but have larger interpolation errors; this is because 
(tg - t 2 ) is larger than the time intervals used for linear interpolation in 
the regular schedule, and Gibb's method is more accurate than linear 
interpolation. 


TABLE 3.- ALTERNATE SCHEDULES FOR TYPE 5 USING ONLY 
LINEAR INTERPOLATION 


Time 


^2 

^3 

t4 

^5 

^6 

ty 

Angle measured, case 1 

- 

s 

Tl 

T 2 

Tl 

S 

- 

Angle measured, case 2 

_ 

Yl 

T2 

S 

T 2 

1 

- 
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Computation of Velocity 


It is a well-known principle of celestial mechanics that if the position 
vector on an orbit is known at two different times, the velocities can be cal- 
culated. In the two-body case this determination of velocity is known as 
Lambert's problem, and a number of methods of solution are in the literature. 
However, these computational procedures are valid for large time intervals and 
are fairly complex. For shorter time intervals Gibb's method (appendix B) 
will provide a reasonably accurate solution with much less computation and 
allow the inclusion of perturbing gravitational forces. 

Two separate sets of measurements of any of the five types already 
discussed will provide two independent position fixes. With types 4 and 5 it 
is also possible to obtain two position fixes by extending a single measure- 
ment schedule to include one additional measurement each of a star- limb angle 
and the subtense angle. The use of the two independent position fixes and the 
extended schedules is discussed in the following paragraphs. 

Two independent position fixes- If r(tx) and r(t 3 ) have been determined 
by two sets of any of the five measurement types, then the velocity v(t 2 ) at 
some intermediate time t 2 can be computed by Gibb's method. It is shown in 
appendix B that if it is assumed that t 2 = 0, then 


v(t2) = 


t3 " ti 


r(t 3 ) - r(ti) - 


A(ti) 


(2t; 


2tit3 - ti2) 


ACtg) 

6 


(t3^ 


2tit3 - 2ti^) 


( 10 ) 


while r(t 2 ) is given by equation (7). These formulas can also be used if the 
epoch t 2 is chosen to coincide with the beginning or end of the interpola- 
tion interval by letting t 2 = tj = 0 or t 2 = t 3 = 0. If, for example, the 
second choice is used equation (7) reduces to 


r(t2) = r(t3) 


and equation (10) becomes 


v(t2) 


r(ti)-r(tq) ti - 

[A(ti) + 2A(t3)] 


( 11 ) 


The choice of epoch requires less computation but is less accurate, as will 
be shown in the section on Error Analysis Methods. 
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Extended schedules- The schedules of types 4 and 5 can be extended to 
provide enough information for computing both r and v from a single set of 
measurements. The measurement schedules and the necessary computations for 
type 4 are shown in table 4. 

TABLE 4.- EXTENDED SCHEDULE FOR TYPE 4 MEASUREMENTS 


Time 


t2 

t3 

t4 

^5 

Angles measured 
Computed from 
measurements 
From interpolation 
(eq. C8)) 

Computed from measured 
and interpolated 
quantities 
From interpolation 
(eqs. (7) and (10]) 

Y ^ and S 
rCti) , 
'('1 (ti) 

Y2 arid S 

r(t2) > 

cos l|Ji Ct 2 ) 
r(t2) 

r(t3) ,v(t3) 

Yi and S 
r(t4) , 

'I'l (tq) 
cos i|i2('t4) 

r(t4) 

Y 2 and S 
r(ts) , 

’I^2(t5) 


Note that this is essentially the same procedure as is used for two 
separate position fixes, the difference being that the extended schedule uses 
some of the measured angles twice. Reusing the data makes no difference in 
the position fix computations, but, as will be seen in the next section of the 
report, it does affect the error analysis. 

The extended schedule and associated computations for type 5 measurements 
are shown in table 5. In this case only 9 measurements are required compared 

TABLE 5.- EXTENDED SCHEDULES FOR TiPE 5 MEASUREMENTS 


Time 

tl 

^2 

I 3 

t4 

15 

te 
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Angle measured 
Computed from 

S(ti) 

Yl (t2) 

S(t3) 

Y2(t4) 

S(ts) 

Yl (te) 

S(ty) 

Y2(l8) 

S(tg) 

measured angle 
From linear 

r(ti) 


r(ts) 


r(ts) 


^(ty) 


r(tg) 

interpolation 
(eq. (9)} 


r(t2) 


r(t4) 


r(te) 


r(t8) 


From measured and 










interpolated 

quantities 


(I 2 ) 


1^)2 ( 14 ) 


4 >i (te) 


'Pz (ts) 


From interpolation 
(eq. (8)) 

From measured and 




>l'l (t4) 


'P2 (le) 




interpolated 










quantities 
From interpolation 




r(t4) 


r(tG) 




(eqs. (7) and 





r(ts) 





(10)) 





v(ts) 
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to the 14 needed for two independent position fixes, but the amount of 
computation is the same. As in the type 4 measurements the error analysis is 
substantially different from that for the regular schedule. 

This approach of making only one set of observations to determine both 
position and velocity is discussed further in the section on error analysis 
and in appendix C. However, as will be shown later, the procedure appears to 
have limited practical value because of accuracy limitations. 


ERROR ANALYSIS METHODS 


The previous section presented equations for finding the position and 
velocity vectors from appropriate sets of measured angles. This section deals 
with the equations used to evaluate the errors produced in the computed posi- 
tion and velocity by measurement errors and by the use of linear interpolation. 
The errors in position due to measurement errors are discussed first and then 
those in velocity. These errors are evaluated to a first-order approximation 
by differentiating the position fix equations from the previous section with 
respect to the measured angles. The errors due to the use of linear interpo- 
lation are then dealt with briefly. The equations used for evaluating the 
measurement errors are derived in appendix C and those for the interpolation 
errors in appendix D. 

The measurement error in each angle is assumed to consist of a bias and 
a random component. It is assumed that the random component is gaussian with 
zero mean and independent of the random components of errors in the other 
measured angles, but errors in Yi and Y 2 have the same standard deviations. 
The resulting errors in position and velocity are evaluated on the basis of 
the standard deviations of their components. The bias errors are discussed in 
detail in the following paragraphs. 

Data presented in reference 5 indicate that even after careful calibra- 
tion there may be fairly large mean errors in a series of measurements taken 
over a short period of time. It is customary to treat such unknown biases as 
random variables in order to determine the rms error for the ensemble of mea- 
surements using different instruments and operators at different times. In 
this study, however, it was desired to assess the relative importance of the 
bias errors in a given practical case where one instrument is used by the same 
operator over a short period of time. Therefore the bias in each measured 
angle is assumed to be a constant, and the resulting errors in position are 
determined directly. 

The biases are separated into two parts, one part due to irradiance or 
other uncertainty in planet diameter and the other part due to instrument 
inaccuracy, including operator error. This separation is made because it is 
pointed out in reference 5 that irradiance is a major source of bias error 
which can be corrected for only approximately at the present time.^ It can be 

^Another important source of bias error is the window- induced error men- 
tioned in reference 5. This error can be corrected for by the method 
described in reference 10 and will not be considered here. 


13 



seen from sketch (c) (taken from ref. 5) that when y is measured to the near 
limb, it is reduced by the irradiance bias bj. On the other hand, if y is 
measured to the far limb, it is increased by bj. 



Sketch (c) 


Data presented in references 5 and 7 indicate that a zero calibration 
error may account for a fairly large part of the instrument bias. This cali- 
bration error, which would be a constant for a given instrument line of sight, 
will be discussed with the error equations. For numerical evaluation of the 
effects of the bias errors, to be discussed in the section on results, the 
signs of the individual instrument biases were chosen to give the largest 
possible error. 

The treatment of these bias errors is pessimistic, and quite simplified, 
but it does show how such biases affect the accuracy of the trajectory 
determination and indicate the need for better calibration. 

The equations used for evaluating the position errors produced by 
measurement errors are presented next. 


Measurement Errors in Position 


The measurement errors are evaulated to a first-order approximation by 
differentiating the position fix equations with respect to the measured angles. 
The differentials of the angles are replaced with the various error components 
described earlier as follows: 


dRA = + bpA 

dB = + bg 

dS = 6g + bg + bj 

dy^ = 6i + bi - kibj 

dy^ = 62 + b2 - k2bj 


( 12 ) 
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The ki are defined by equation (1) while 6 and b subscripted by the 
measured angle are the random errors and instrument bias in that angle. The 
irradiance bias is denoted by bx- 

If the spacecraft-planet vector r given by equation (2) is 
differentiated, then 


/uA 

dr = dr I U2 ) + r I du2 j 
\u3/ \du3/ 


Since u and du _are orthogonal dr can be separated into two components, dr 
along r and r du normal to it. The magnitude of r du can be written 


r|du| = r/dux^ + du 2 ^ + du 3 ^ (13) 


The position errors due to measurement biases are found by setting the 
random errors in equations (12) equal to zero. The radial bias error b^. is 
obtained by substituting the resulting differentials into the appropriate 
expression for dr while the normal bias error bj^ is obtained from 
equation (13) . 

The random errors are evaluated on the basis of the standard deviations 
of the components of dr. These standard deviations, denoted for the 

radial component and for the normal component, are given by 


= /e (dr^) - br^ 

= >4‘^E(dui^ + du2^ + du3^) - bjx^ 


(14) 


Equations (13) and (14) show that 


E|dr| = /(a ^2 + ^^ 2 ) + (bj.2 + b„^) 

and the magnitude of dr, which is the vector error in r, is the same as that 
of a variable having a constant bias of magnit ude /b^.^ + and a random 

component with standard deviation /a^2 + a^2. The following paragraphs give 
the equations, derived in appendix C, for br> bj^, o^., and for the differ- 
ent types of measurements. The standard deviation of each of the random 
errors in equations (12) is denoted by where the subscripts are 

the same except that it is assumed that oj = 02 = 
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Type 1 measicpements- The biases are 


bj. = - (bj + bs)r/r2 - 1 


b^ = rv43jj2 + cos^ D 


CTr = cTci‘>'r2 - 1 




Type 2 measurements- The bias errors are 


bj. = -(bj + b5)r/r2 - 1 


bn = [(bi + kibs)^ + (ba + k2bg)2 - 2 (bi + kjbs) (b2 + k2bs)cos 


It was stated earlier that the signs of the biases were chosen to give the 
largest possible error; therefore if it is also assumed that 
|bil = |b 2 | = |bs| = b, then 


|bj.| < (bj + b)r/r2 - 1 
|bnl + Icos B|) 


The standard deviations of the random errors are given by 


agr/r2 - 1 


■ ^ n \2aJ^ + 2aQ^ (1 + kik 2 cos B) 
sin BY * ^ ^ 


For this type of measurement there are singularities in both b^ and when 
B = 0 or 180°. Thus, poor accuracy would be expected near these values of B. 

Type 3 measurements- In this case the errors are the result of errors in 
Y 2 j Y 2 j b, and 
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r / r ^ - 1 


bj. = 2~" ~ [^1^1 Cl + Pl^ k2b2(l - P^) + 2bgP2 - 2bi] 


sin B 


I »ii>i 


^ 2 ^ 2 ^^ 


[(1 + Pj2) + kik2(l - Pi^)cos B] 


+2bB(kibi - k2b2)P2p2fl " cos B) + 2bB^p2^(l - kik 2 cos B) 


1/2 


If |bi| = |b 2 | = b and the signs of the biases are chosen to give the largest 
possible error, then 


bj.| < r/r2 - 1 [bi + b (1 + Ip^I) + |bgp2|] 


l*n 




sin B 


{b2[(i + + kik2(l - Pj^)cos B] 


+4b|bg| |p^P 2 Cl - kik 2 cos B) | + 2bB^P2^(l " kik 2 cos B) } 


1 /2 


( 20 ) 


The standard deviations of the random errors are 


- 1 (1 + p^^) + 


CB ^ P 2 ^ 


O-n = 


n " Sin B 


{a 2[(i + p^2) + kik 2 (l - Pi2)cos B] + 2aB2p2^(l - kik 2 cos B)}l'2j 


In equations (20) and (21) p^ and P 2 are defined by 


( 21 ) 


sin(kn|;i - k2<{)23 (1 ^^ 1^2 cos B) 

~ sin(kx4)i + k2<{'2) (1 ” ^^1^2 cos B) 


and 


sin i|)x sin i /»2 ® 


2 sin(k 2 'f'i + k 2 'j' 2 ) (1 “ ^^ 1^2 cos B) 


If ki = k 2 , it can be shown that 
1 + px^ + kxk2(l - px^)cos B 


sin^(ij;x - 1P2) 


sin^ B 


1 - cos B sin^(i|;x + ^2) (1 " cos B)' 
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and 


P 2 ^(l - kj^k 2 cos B) sin^ sin^ ip2 

sin^ B + i|^2) (1 “ B) 

In that case for type 3 has singularities for B = 0 and Jpi + ^2 ~ 180°, 

but none for B = 180°. Similarly, if kj = -k 2 there are singularities for 

B = 180° and = ip 2 but none for B = 0. 

Type 4 measurements- The equations for br and bj^ for type 4 are 
identical with those for type 2 and are given by equations (18). Also, since 
r is determined from a single measurement, Oj. is the same as for type 2, and 


Oj. = - 1 

"" s ' i n B |f '^“1“ ^/r^ - 1 ctn (2k2 cos Bi + kj] 

( 22 ) 

The equation for was derived under the assumption that is the angle 

interpolated. If ijj 2 is interpolated instead, it replaces in the 

equation and k]^ and k 2 are interchanged. 

In addition to the singularities for B = 0 or 180° noted for type 2, 
Gibb's method produces a singularity when the interpolated angle ^2 

zero. This singularity is in the measurement error and should not be confused 
with the interpolation errors discussed later in the report. 

lype 5 measurements- The bias errors for this type of measurement are also 
given by equations (17), but is was pointed out earlier that a large part of 
the bias errors appear to be due to zero calibration errors. If the biases in 
Y and S are attributed entirely to this source, then equations (17) for 
type 2 biases remain the same. However, since all angles for type 5 are 
measured with the same instrimient b^ = b 2 = bg = b and 



|br| < (bi + b)r/r2 - 1 

- 2(1 + ki + k2 + kik2)cos B]^'^ 


(23) 


These equations were used for type 5 in order to assess the effects of the 
zero calibration errors. Note that if k 2 = k 2 = 1, there is no singularity 
for B = 0, and if ki=k 2 = -l, then b^^ = 0. 

The standard deviations of the random errors are 
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( 24 ) 


Type 5 has singularities at the same points as type 4, that is, at B = 0, 
B = 180°, and \pi = 0. Since types 4 and 5 both have terms containing ctn ipx 
where \pj_ is the angld interpolated, it will be assumed that interpolation 
will always be carried out with the larger of the two angles. 

Note, also, that Oj. for type 5 is reduced by compared to that for 
types 2 and 4 because two measurements of S are used in determining the 
range. 

The two alternate schedules for type 5 (discussed in the section on 
position fix equations) which use only linear interpolation have nearly the 
same position errors due to measurement inaccuracies as type 2. These 
schedules can be evaluated on the basis of data from type 2. 


Extended schedules- It was pointed out in the section on the computation 
of velocity that the velocity could be determined by adding one measurement 
time to the type 4 schedule or two to the type 5 schedule. The resulting bias 
errors are unchanged, but the standard deviations of the position errors 
become, for type 4, 



and for type 5 
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/sCr^ - 1) 
ar = asr J ^ 

" 2 sin B ’I'l + 2 ctn^ \p 2 


+ cos ctn \p2 cos B] + s/r^ - 1 (ki ctn ^2 ■*■ ^2 ctn i|^i)cos B 

1 1/2 

- 15kik2 cos B + 20 } 


( 26 ) 


Note that for both measurement types, is reduced because more measurements 
of S are used in determining r. On the other hand, singularities are 
present for both ijji = 0 and 1 P 2 = 0. 


Velocity Errors 

It is shown in appendix C that, for most conditions, the third-order 
terms in equation (10) may be neglected in computing the measurement errors in 
v(t 2 )’ That is, we can assume 


dv(t2) 


drCtg) - dr(t;^) 
^3 " 


( 27 ) 


Since the biases are considered to be constant over the observation period, 
the portions of df(tj^) and df(t 3 ) due to bias errors are equal. Therefore, 
to a first approximation dv(t 2 ) is due only to random errors, and the 
standard deviation, Oy of |dv(t 2 )| is 


/ 


2 (Oj-^ + Oji^) 

^3 - tl 


/ 


^fi At 


( 28 ) 


Thus, to a first approximation the standard deviation of random error in the 
velocity is inversely proportional to the time between the two position fixes. 
The standard deviation of error in r(t 2 ) (see appendix C) is 


Cp(t2) 


/(Or^ + On^)(tj^ + 

T 3 - tl 


( 29 ) 


This error is minimized if -ti = t 3 = At, in which case 


Cp (t2) 
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and 


(to) 

^v(t2) = - (30) 

If Gibb's method is used to obtain v(t 3 ) instead of v(t?) , then r (ts) is 
assumed to be the correct position, and Op increases to since 

the information obtained at t^ is not used. However, Oy is the same as 
when V is computed at t 2 - Equation (30) is roughly the same (usually 
within a factor of two) for the extended schedule versions of types 4 and 5. 


Linear Interpolation Errors 

It was pointed out in the description of type 4 and type 5 measurements 
that it is necessary to use some linear interpolation with type 5. If one 
component of position (in other words r and either or i|^2) is known at 
two times a more accurate interpolation, such as Gibb's method, may be used. 
Since linear interpolation is simpler than Gibb's method, and in some cases is 
necessary, we would like to know when the resulting errors are small enough 
compared to the measurement errors to be acceptable.^ 

The bias and random errors are a function only of the geometry, but 
errors due to linear interpolation of r or y depend on the orbit as well. 

For this reason only upper bounds of these errors were considered. The inter- 
polation error for r was evaluated by taking the maximum magnitude of the 
second-order term in a Taylor series expansion of r. This error increases as 
r decreases, and the minimum value of r was foun d for wh ich the 

interpolation error is less than some fraction, M, of Ogr/r2 - 1 (i.e., of 
a-j. for type 4). In order to ensure that the second-order term is a valid 
measure of the interpolation error, an expression was also found for the mini- 
mum value of r such that the third-order term is negligible compared to the 
second-order term. 

A conceptually similar, but less direct, approach was used for the errors 
due to linear interpolation of y. The minimum value of r was found for 
which the interpolation error is less than M times /a^2 + for type 5. 

Tne error due to interpolation of y is in the normal direction, but the 
radial component is often the dominant part of the measurement error and 
therefore should be included in determining the relative importance of the 
interpolation error. The value of r^in in this case is a function of B as 
well as r. 

The values of rmin for which the interpolation error is less than M 
times the appropriate standard deviation are presented in the section on 
results. Appendix D contains the derivations of the equations used and the 


^There is also an interpolation error due to the use of Gibb's method 
which is much smaller than that due to linear interpolation. This error, 
which cannot be analyzed by linear perturbation methods is discussed briefly 
in the section on results . 
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values of rmin for neglecting the various higher order terms. The latter 
results indicate that in most practical cases the second-order terms are a 
valid measure of the interpolation error. 


RESULTS AND DISCUSSION 


The standard deviation of position error, and the total bias 

error, were computed for all five types of measurements with 

various values of RA, D, and <J> at several different ranges from the planet. 
At each value of r sufficient cases were computed to cover the range of 
angles considered (-90° RA ^ 90°, 0 ^ D ^ 90° and 0 ^ <j) ^ 90°) . 


It was found that, except for type 1, the errors depend strongly on the 
location of the projection of r onto the X-Y plane relative to the unit 
vectors to the two stars. Three different star-planet configurations, two of 
them extreme possibilities for location of the projection and the third an 
intermediate value, are illustrated in figure 2. In the first configuration 



(a) First configuration, A = 0/2. 




Figure 2.- Star-planet configurations. 

(fig. 2(a)) the projection of r on the X-Y plane lies midway between the 
two stars. In figure 2(b) the projection of f is in the direction of one of 
the stars. In this case, star number 1 and RA = 0, but it could be in the 
direction star number 2 instead. In figure 2(c) the projection of f lies 
well outside the smaller angle between the two stars. 

The values of + Oj^2 and + bj,^ obtained for these three 

geometric configurations are discussed separately in the next two sections 
after which the errors caused by linear interpolation are compared with those 
due to measurement errors. 


Random Errors 

The standard deviations for the random position errors were computed 
assuming the same standard deviations of error in each angle measured, that is. 
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Ora = Od ~ °y ~ ~ ~ except for a special set of data for type 3 

measurements. Tests made on the trisextant at Ames Research Center indicated 
that it may not be possible to measure B as accurately as the other angles. 
Therefore, an additional set of data was obtained for type 3 measurements with 
Or = 10c. 

The value of + (j^2 £qj. type 1 measurements is presented in 

figure 3. The standard deviation in units of (planet radii)/ (arc sec) is 
plotted as a function of the declination D for four ranges. The correspond- 
ing values of or (dashed lines) are also presented. Note that as r 
increases the importance of Cj- in the total standard de viation bec omes 
greater, and for r ^ it accounts for nearly all of /oj.2 + , Since 

Oj. is independent of D, there is negligible variation of the total standard 
deviation with D at the larger ranges. It can be seen from equations (16) 
that these results are to be expected and that for large values of r the 
type 1 error may be regarded as consisting only of the radial component. 

IO-3|- 


r=IO planet radii 


1 10 


-4 


I0‘5 


<^r 




r=2 

r=2 


-r=yr 

-r=yr 


I I I I I I I 1 I I 

0 10 20 30 40 50 60 TO 80 90 

Declination, deg 


Figure 3.- Standard deviations of random error for type 1 
measurements. 

The errors from measurements of type 2, 3, 4, or 5 are functions not only 
of the range but also of the relative location of the stars and planet and 
whether the star- limb angles are measured to the limb near the star or the far 
limb. The ratios of the standard deviations for these types of measurements 
to those for type 1 were computed in order to facilitate comparison between 
the various methods. The effects of the relative locations of stars and 
planet were evaluated by using the three different configurations (fig. 2) dis- 
cussed earlier. Two different values of r, /To, and 10, were found to be 
sufficient to show the effect of changing range. The results of measuring yi 
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and Y 2 near and far limbs were accounted for by computing the 

22 cases described in table 6. 

TABLE 6.- CASES FOR WHICH ERROR RATIOS WERE COMPUTED 


Case 

Measurement type 

Description 

2a 

2 

and y 2 to same limb 

2b 


Yi and Y 2 to opposite limbs 

3a 

3 

Y 2 and Y2 to same limb, Og = 0 

3b 


Yi and Y 2 to same limb, Og = 10 a 

3c 


Yi and Y2 to opposite limbs, og = a 

3d 


Yi and Y2 to opposite limbs, Og = 10 

4a 

4 

Yi and Y2 to near limb 

4b 


Yi and Y2 to far limb 

4c 


Yi to near limb, Y2 to far limb 

4d 


Yi to far limb, y 2 to near limb 

4e-4h 

4 

Same as 4a- 4d with extended schedule 

5a-5h 

5 

Same as 4a- 4h 


a 


The resulting ratios of types 2, 3, 4, and 5 to those for 

type 1 were computed for each of the geometric configurations and each of the 
two values of r. These results are presented graphically in figures 4 
through 12. The curves for type 2 are presented with each of the other types 
separately for ease of comparison. Only the maximum and minimum values from 
cases 4(a) through (d) are shown and are labeled as type 4, while the maximum 
and minimum values for the extended schedule are labeled 4X. An equivalent 
set of data is presented for type 5 measurements, and the curves are labeled 
5 and 5X. 


The data for the configuration with RA = cj)/2 (fig. 2(a)) ar e p resented 
in figures 4, 5, and 6. Three values of (j) are used with r = /To, and one 
is also used with r = 10 (figs. 4(c), 5(c), and 6(c)) in order to show the 
effects of range. 

Figure 4 shows the ratio for types 2 and 3. For this con figuration 
B = 180° when D = 0 and decreases as D increases so that /Oj.2 + for 

type 2 is infinite when D = 0, reaches a minimum at the value of D for 
which B = 90° and then increases again. This effect is most noticeable for 
the smallest value of <p since the minimum value of B increases with (f) . 

On the other hand, the curves for types 3a and 3b are at a minimum for D = 0 
and increase rapidly for larger values of D. When RA = <j>/2, we can see that 
'('1 = causing (as was pointed out in the section on Error Analysis Methods) 
the errors for types 3c and 3d to be infinite, and those curves have been 
omitted. For (p = 30° and D ^ 5°, type 2 errors are less than twice those 
for type 1, and the ratio is smaller for larger values of cp. Also, compari- 
son of figures 4(b) and 4(c) shows that as r increases, the ratios for both 
types decrease and the spread between the curves of the same type becomes 
smaller. From the equations for the standard deviations we would expect the 
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Declination, deg 

(a)0 = 3O°,RA= 15°,r = v^ 


O 10 20 30 40 50 60 70 80 90 

Declinotion, deg 

(b) = 60° , RA = 30° , r = v/To" 
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Declination, deg 
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0 10 20 30 40 50 60 70 80 90 

Declination, deg 


(c)0 = 6O°,RA = 3O°,r=lO 


(d) 0 = 90°, RA = 45°, r = v/To 


Figure 4.- Ratio of \/o^ + for types 2 and 3 with RA-0/2 to for type 1. 


sensitivity to changes in r to be greater for small values of r. Therefore, 
for (p ^ 30°, r > v'To^ and any given value of D. One or both of these 
measurement types (type 2 and type 3) will provide a standard deviation less 
than twice that for type 1. 



The ratio for type 4 (fig. 5) follows the same trends for changes in 
r or (j) as the corresponding data for type 2. However, type 4 is generally 
poorer than type 2, particularly for small values of D. The use of the 
extended schedule improves the results somewhat and produces lower ratios (to 
type 1 values) than type 2 for the larger values of D. This improvement is 
due mainly to the increase in the number of measurements used. 




\ 

\ 
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(a)(/) = 30°,RA= 15°,r = v/T0 


(b) 0 = 60° , RA = 30° , r = vTo 



r i I I 1- fill II 
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Declination, deg 

(c)0 = 6O°,RA = 3O°,r= 10 
Figure 5.- Ratio of for type 2 and 



(d)0 = 9O°,RA = 45°,r-v/To 
4 with RA = 0/2 to for type 1 . 
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The results for type 5 are given in figure 6, and the trends with changes 
in r and ip are the same as those for types 2 and 4. Type 5 has substan- 
tially smaller errors than type 4 and less spread between maximum and minimum. 
The singularities at D = 0 have a greater effect on type 5 than on type 2, 
but for D > 25° the type 5 errors are smaller than those for type 2. This 
improvement in performance is probably due to the increased number of obser- 
vations for type 5, and for the same reason the use of the extended schedule 
urovides a small improvement over the regular schedule for type 5. » 



(a)b = 30°,RA= 15-,r=VrO (b) b = 60°, RA = 30°, r = VTO 


6 - 
5 - 

4 - 



0 10 20 30 40 50 60 70 80 90 

Declination, deg 



(c)0 = 6O°,RA = 3O°,r= 10 


(d) 0 = 90°, RA = 45° , r = VTo 


Figure 6.- Ratio of Vo^TTo^ for type 2 and type 5 with RA = 0/2 to for type 1. 
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The corresponding data for the configuration (fig. 2(b)) where the 
projection of f on the X-Y plane is alined with one of the stars are 
presented in figures 7, 8, and 9. The errors for this configuration show the 
same trends with r and 4> as shown in figures 4, 5, and 6; namely as r or (j) 
increases j the ratios of the standard deviations to those for type 1 become 
smaller, and the spread between the maximum and minimum curves is reduced. 

The effect on types 2 and 3 of changing the geometric configuration can be 
seen’ by comparing figure 4(a) with figure 7(b) and figure 4(b) with fig- 
ure 7(d). These pairs of figures have the same values of <}> and r and the 
curves for type 2 are nearly same except for small values of D. The type 2 



Declination , deg Declination, deg 


(a)d= 15°,RA = 0, r = VIO (b) d = 30°, RA = 0, r = VTO 



Declination, deg Declination, deg 

(c)0 = 3O°,RA = O,r= 10 (d) 0 = 60°, RA = 0, r = Vib 

Figure 7.- Ratio of \/apTa^ for type 2 and type 3 with RA = 0 to + On^ for type 1 . 
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errors for RA = (j>/2 increase quite rapidly for D < 20° while there is 
little noticeable increase for RA = 0. On the other hand, the errors for 
type 3 are much larger for RA = 0 than for RA = but the combination‘s 

of types 2 and 3 will still provide a standard deviation less than twice that 
of type 1 for all values of D if cj) ^ 30° and r ^ /To”. 

When the regular schedule is used, the errors for type 4, shown in 
figure 8, are not affected as strongly by the singularities near D = 0 as 



(a)<l>= 15°,RA = 0, r = VT0 



gl 1 I I I 1 I 1 I I 
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(c)0 = 3O°,RA = O,r= 10 
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(d)(^) = 60°,RA = O,r = \/T0 


Figure 8.- Ratio of Va^^~+a^ for type 2 and type 4 with RA = 0 to for type 1 . 

*SNote that if RA is exactly zero then lim B = 90°, and a-p for type 2 

D 0 

(eq. (19)) is bounded. However, for RA slightly different from zero, Oj- 
will approach infinity and it would be necessary to use type 3 for very small 
values of D. 
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with the configuration where RA = <j)/2. On the other hand, the singularities 
have greater effect on the results for the extended schedule when RA = 0. 
This result is to be expected since , which goes to zero with D is not 

interpolated with the regular schedule, but it must be used with the extended 
schedule thereby introducing ctn ifij , into the equation for Ojj. 


Figure 9 shows the same data for type 5 measurements, and the same 
changes with r and (() are noted as when RA = (j>/2. The change in geometric 
configuration from RA = ifll to RA = 0 has the same general effect on the 
type 5 results as on those for type 4, and for both configurations the ratios 
for type 5 are considerably lower than those for type 4. 
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Figure 9.- Ratio of sja^ + Ojx for type 2 and type 5 with RA - 0 to for type 1 • 
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The data for the configuration (fig. 2(c)) where the projection of r 
onto the X-Y plane lies well outside the smaller angle between the two stars 
are presented in figures 10, 11, and 12. The ratios for types 2 and 3, which 
are shown in figure 10, are much larger than those for the other two configu- 
rations. Comparison of figures 10(a), 10(b), and 10(d) with figure 7(a), 
which has the same values of <p and r, shows that the ratios for type 2 are 
the same at large values of D, but the ratios in figure 7(a) are becoming 
relatively much smaller as D decreases. The ratios for type 3 start at 
larger values for D = 0 but do not increase as rapidly as those for the 
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other two configurations. Note that the ratios for types 3a and 3b are off 
the scale of the graphs in all cases, and in figure 12(d) the curve for 
type 3c has been omitted for the same reason. As with the other two configu- 
rations the ratios decrease with increasing r, but configurations of this 
type should be avoided for r < 10. 

The ratios for type 4 with the third configuration are shown in figure 11, 
and are quite large in all cases except at very large values of D. Note that 
for RA = 90° (fig. 11(d)) the error for the regular schedule is independent 
of which limbs and y^ are measured to so that only one curve is shown. 



Dedinotion , deg Declination, deg 

(a)0= 15°,RA = 30°,r = vT0 (b)0 = 15°, RA = 60°, r = VTo 
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(c)0= 15°,RA=60°,r= 10 (d) 0 = 15°, RA = 90°, r = vTF 

Figure 1 1 Ratio of for type 2 and type 4 with RA > 0 to for type 1 . 
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Figure 12 presents the same data for type 5 measurements, and, again the 
ratios are considerably larger than for the other two configurations. However, 
for r > 10 and D > 25° the type 5 errors are comparable to those for type 1. 


The preceding results show that measurements of type 2, 3, or 5 can have 
accuracies comparable to that of type 1, provided it is possible to choose the 
stars properly. Type 4 measurements appear to be decidedly inferior to the 
others and will not be considered further. For type 2 the stars should be 
chosen so that B is as near 90° as possible, while for type 5 (regular 
schedule) either ipi or i {)2 as well as B should be as near 90° as possible. 
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These conditions are satisfied best when the projection of r onto the X-Y 
plane lies near one of the stars. For type 3 the best results are obtained if 
the two stars and the planet are nearly colinear with the projection of f 
onto the X-Y plane midway between the two stars. In all cases the results 
are poorer when the projection of r lies well outside the smaller angle 
between the two stars. 

The practicality of making a suitable choice of stars has not been 
investigated in detail here, but an examination of the star chart in refer- 
ence 6 indicates little difficulty if a 90° field of view is available. Reduc- 
tion of the field of view to 55° makes it more difficult to locate suitable 
stars, and in the case of type 5 increases the minimum obtainable error since 

If we ass ume that the proper choice of stars can be made, then the ratio 
of /a^2 + Oj^2. for types 2, 3, and 5 to that for type 1 can be held to less 
than 2. As r increases, this ratio approaches unity over an increasing 
range of D because the radial component becomes the dominant part of the 
error for all types. In view of these facts it seems reasonable to assume 
that the total random error for any of these observation types can be approxi- 
mated by the radial error from a single subtense angle measurement. There is 
some improvement in the results for type 5 with the extended schedule but the 
change is not enough to invalidate this rough approximation. 

This approximation used with equation (16) for r = 200 indicates position 
errors with a standard deviation of about 0.2 planet radii for each arc second 
standard deviation of instrument error, or, in other words, an error of 1 or 
2 percent for a 10 arc sec instrument. Thus, within the present state of the 
art, the measurements considered in this study are useful to a maximum range 
of about 200 planet radii from the central body. 


Bias Errors 

The computation of the bias errors is complicated by the fact that 
different biases may either add or subtract. Therefore, the upper bound of 
/bj.2 + was computed by assuming signs on the different biases which gave 

the largest total error. The instrument biases in y, S, A, and D were 
assumed to have the same magnitude b, while the bias in B and that due to 
irradiance were assumed to be either b or 10b. It was assumed that b is 
equal to the standard deviation a of the instrument errors discussed in the 
previous section. 

The ratio of the biases for type 1 measurements to the la random error 
for type 1 were computed and are plotted in figure 13. These ratios are 
nearly constant and are approximately (1 + lbi|/|b|). This approximation 
becomes more accurate as ]bil or r is increased. 

The ratios of /bj.2 + b^2 for the remaining measurement types to Cj. for 
type 1 were computed for the same cases as the random error and are plotted in 
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figures 14 through 16. Figure 3 shows 
that, for the values of r considered, 
the type 1 random error is almost 
entirely in the radial direction, and 
it is reasonable to compare these 
ratios directly with those in fig- 
ures 4 through 12 . Note that , as was 
pointed out earlier, the type 5 biases 
were computed under the assumption 
that the entire bias is due entirely 
to zero calibration error; actually 
the ratios for types 4 and 5 will lie 
somewhere between the type 2 and 
type 5 curves, depending on the rela- 
tive magnitude of the zero calibration 
error. Except for one case (fig. 16(d)), 
the ratios for bj = 10b are about 
five times those for bj = b and the 
value of bj is not indicated on the 
figures . 

The ratios for the configuration 
with RA = 4>/2 (fig. 2(a)) are pre- 
sented in figure 14. The curves show 

the same general characteristics as those for the random errors in figures 4, 

5, and 6; that is, types 2 and 5 have singularities at D = 0 while type 3 
has a singularity at D = 90°, and the ratios for types 3c and 3d are infinite 
because of the singularity when ipi = tp 2 - Increasing r from /To^ to 10 
(figs. 14(b) and 14(c)) reduces the ratios for bj = 1 but has little effect 
for bj = 10. Note that the smallest value of /b^2 + bj^2 , type 5b with 
bj = 1, is twice as large as for type 1, while if bj is increased to 10 

the ratio increases to 11. The biases due to zero calibration error are 
reduced by measuring one or both of the star- limb angles to the far limb, but 
the reduction is significant only at small values of D. 

Figure 15 presents the corresponding data for the configuration where 
RA = 0 (fig. 2(b)), and, again, the curves follow the same general pattern as 
the corresponding ones for the random errors (figs. 7, 8, and 9). In fig- 
ure 15(a) the ratio for type 2 is much larger than that for types 5c and 5d 
which is, in turn, much larger than the other ratios for type 5. This fact 
serves to show that the bias errors increase in magnitude and spread between 
various cases in the same manner as the random errors as a singularity is 
approached; in this case B becomes small with increasing D. 

The configuration with the projection of r well outside the smaller 
angle between the two stars (fig. 2(c)) produced the ratios illustrated in 



Figure 13.- Ratio of bias errors to la random errors for 
type 1 measurements. 
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Figure 15.- Ratio of bias errors for types 2, 3, and 5 with RA = 0 to Oj- for type 1. 
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figure 16. As in the case of the random error (figs. 10, 11, and 12) the 
biases for type 2 and types 5c and 5d are much larger at lower values of D 
than those for (j) = 15° and RA = 0 (fig. 15(a)). The other biases, types 3, 
5a, and 5b, are generally a little smaller than for the other configurations. 

The upper bounds on /bj.2 + bj^2 presented in these figures are only rough 
approximations and represent worst possible cases. However, they do show that 
in many instances a bias as large as the standard deviation of the random 
error can be the dominating source of position error. Even if we attribute 
the entire bias to zero calibration error and use type 5b measurements so that 
bn = 0, the bias error is still about twice the random error. An increase in 
bj causes a large increase in the bias error and every effort should be made 
to correct for irradiance bias; even a crude correction would be helpful. 

As in the cases of the random errors, the bias errors for measurement 
types 2, 3, and 5 can be reduced to about the same value as those for type 1 
by the proper choice of stars. Furthermore, the same general conclusions as 
to which stars should be chosen to minimize random errors also apply in the 
case of the bias errors. 

In order to draw any further conclusions with regard to the relative 
importance of bias and random errors, it is necessary to make additional 
assumptions regarding the numerical values of the biases and standard 
deviations. For example, the instrument bias error from the Gemini flight 
described in reference 8 averages about 2.25 arc sec or about one third the 
average standard deviation. The irradiance error for the Moon as seen from 
the Earth is about five times the instrument bias error (ref. 5). 

We can use these numbers with equations (15) and (16) to show that the 
position bias errors will be approximately twice the standard deviation of the 
random errors, and the total rms error is about 2.2 times the standard devia- 
tion. On the other hand, if the irradiance error could be reduced by calibra- 
tion to about the same value as the instrument bias error, the total rms error 
would be about 1.3 times the standard deviation of the random part. 


Position Errors Due to Linear Interpolation 

The errors caused by the linear interpolation of r and y are dealt with 
separately in this section. The derivations of the equations used and some 
data pertinent to neglecting higher order terms are given in appendix D. 

Interpolation of r- Appendix D gives an expression for rmin? the 
minimum value of r for which the error in the linear interpolation of r is 
less than Ma^., where is evaluated for type 1 measurements. The value of 
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r^iin is a fianction of the astrodynamic constant y, the orbital eccentricity, 
e, the standard deviation of random measurement errors, a, and At, which is 
half the time between the observations being interpolated. The computation of 
rjnin requires the assignment of values to the various parameters involved. 

The value of y for the Earth was used since it gives the largest value of 
rmin major body in the solar system (see appendix C) , 10 arc sec was 

used for o and values of e were used which cover the range of values 
likely to be encountered in manned space flight. 

The data in reference 8 indicate that sextant (single angle) measurements 
of the same angle can be repeated on the average at intervals of about one min- 
ute. Since the alternate measurement of different angles or the simultaneous 
measurement of two angles will require more time, rjj^qj^ was evaluated for At 
of 1, 2, 3, 4, and 5 minutes. 


Figure 17 shows the resulting values of rmin ^ ~ 0.1, plotted as 

a function of At for various values of e. It was found (eq. (D9)) that 

rmin is a function of eAt^/Mo so that the data can be used for other values 

of e, M, a, or At. For example, the 
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is the same as for 
M = 1.0. 


: 1.0 and M = 0.1 
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If At is 5 minutes, M = 0.1 
and e = 0.98, the equivalent of a 
lunar trajectory, then rmin Earth 
radii, while if At can be reduced 


to 1 minute r 


min 


5 Earth radii. 


= .001 


At, mm 


Figure 17.- Minimum radius for which linear interpolation 
error of radius is less than 0.1 a^. 


Although these are rather large dis- 
tances, they are only 4.3 hours and 
1.7 hours, respectively, from perigee 
and r > rmin for most of the tra- 
jectory. This is also true for 
e = 10, in which case rmin ^20 Earth 
radii for At = 5 minutes and 9 Earth 
radii for At = 1 minute. These dis- 
tances are about 1.7 hours and 

40 minutes, respectively, from perigee. Thus, linear interpolation of r will 
produce negligible errors over the major portions of both lunar and interplan- 
etary trajectories. For near-planet orbits, where r 4, interpolation may 
be the major source of error unless e is very small. For example, if 
r = /2, e = 0.1 and At = 5 minutes, the interpolation error would exceed the 
measurement error by about a factor of 50. This factor decreases linearly 
with e and with the square of At, but there are many combinations of r, e, 
and At in the range considered for these parameters for which the inteipola- 
tion is the dominant source of error. This fact does not necessarily preclude 
the use of linear interpolation for near-planet orbits, however, because the 
error due to measurement inaccuracy for the model assumed in the study are 
quite small near the planet. 

InteTpolati-on of y- The position error due to the linear interpolation 
of Y is normal to the radius and contains ctn xjjQ (eq. (D31)), but since the 
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Figure 18.- Minimum radius for whic h error due to 
interpolating 7 is less than . 


radial portion of the measurement 
error is usually dominant, it should 
also be considered in evaluating the 
interpolation error. Therefore, two 
expressions were found (eqs . (D32) 
and (D34)) for which the interpola- 
tion error will not exceed 
M/oj.2 + the value of /oj.^ + a^2 

for type 5 was used since it is 
usually smaller than that for type 4. 

The value of rm-i n in this case 
is a function of sin B as well as 
e, M, a, and At^, and the last three 
terms occur in the expression for 
rmin in the form At^/Ma. Thus a 
single curve represents a family of 
values of M, a, and At, but e and 
sin B enter in a more complicated 
form. 


The resulting values of r^in 
are plotted in figure 18 for two dif- 
ferent values of sin B. For given 
values of e and y, is gener- 

ally larger than the corresponding 
values in figure 17. If the data 
from figure 17 are used in the pre- 
vious example, it is found that for 
the lunar trajectory, with M = 0.1, 
sin B ^ 0.5 and At = 5 minutes, r„,pj^ 
occurs about 7.7 hours from perigee, 
and if At is reduced to 1 minute, 
the time is about 3 hours. For the 
high eccentricity (e = lOj trajectory 
the times are 1 . 8 hours and 45 min- 
utes, respectively. For small values 
of sin B these times become larger, 
but if stars can be chosen to keep 
sin B large, it should be possible 
to use linear interpolation for y 
with lunar and interplanetary trajec- 
tories. For near-planet orbits, the 
values of rjj,in relatively large 

compared to those in figure 17, but 
the absolute errors due to linear 
interpolation at these ranges are 
small and may be acceptable. 

Potential application- The data 
just presented have shown that at 
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sufficiently large ranges the linear- interpolation error is negligible 
compared to that resulting from measurement errors of 10 arc sec standard 
deviation. This fact implies that at such ranges it would be reasonable to 
use the alternate schedules for type 5, which employ only linear interpolation, 
as described in the section on measurement types. It was pointed out in the 
section on error analysis methods that the errors from the alternate schedules 
are about the same as those for type 2, thus comparing favorably with the regu- 
lar type 5 schedule. Furthermore, some error due to linear interpolation of 
r is present even when Gibb's method is used for type 5. The linear- 
interpolation errors would, of course, become relatively more important if the 
instrument errors were reduced significantly or if interpolation intervals 
(2At) were larger than 10 minutes. However, these results make the alternate 
schedule potentially very attractive. 


Velocity Errors 

Once the position vectors have been obtained at two different times, the 
velocity is obtained by interpolation using Gibb's method. Linear interpola- 
tion could also be used to obtain velocity, but, for reasons to be discussed 
later, it is not considered accurate enough. The two major sources of error 
in velocity are the inaccuracy of the interpolation formula and random measure- 
ment errors, the measurement bias errors in velocity being zero. The standard 
deviation, Oy, of the velocity error is approximately (see eq. (28)] 

+ 0^2/ /J At. Since this is the case, the data presented for type 1 posi- 
tion errors in figure 3 can be used for velocity errors by dividing the verti- 
cal scale by /2 At. The ratios of the standard deviations of velocity errors 
for types 2, 3, 4, and 5 to that for type 1 are the same as the ratios for the 
corresponding position errors given in figures 4 through 12. 

The remainder of this section is concerned with the combination of 
velocity errors resulting from random measurement errors and from interpola- 
tion. For this purpose a further simplification, which was pointed out in the 
discussion of random position errors, will be made. It will be assumed that 
/aj.2 + for each position vector, can be approximated by that of the 

radial error arising from a single subtense angle measurement, that is, 

Osr/r2 - i . 

Since the interpolation error increases with At while the random 
velocity error decreases, a value of At can be chosen for which the total 
velocity error is minimized. A rough approximation to this optimum At was 
obtained by trial and error for a few sample Earth orbits which were chosen to 
be circular for ease in calculating the correct r and v. 

Table 7 shows these approximately optimum values of At, along with the 
associated values of Oy, a„, and the errors Av in velocity and Ar in 
position that are due to interpolation. In the table t 2 denotes the center 
of the interpolation interval while the errors at ts apply to either end of 
the interval, Ar(t 3 ) being omitted since it is zero. The error quantities are 
expressed in meters per second and kilometers. The standard deviation of 
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TABLE 7.- COMPARISON OF RANDOM ERRORS WITH INTERPOLATION ERRORS 




Standard deviations 

Interpolation errors from 

r. 

At, 
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min 
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Op(t3) » 
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Ai'(t2) , 

Av(t3) , 

radii 



m/s 

km 

km 

m/s 

km 

m/s 

/2T0 

1 

7.3 

0.31 

0.44 

0.2 

0 

1.3 

/27o 

5 

1.5 

0.31 

0.44 

.4 

6 

25.0 

/IF 

10 

4.8 

2.0 

2.9 

.3 

1 

3.2 

10 

50 

10.3 

22 

31 

.0 

2 

1.8 

Ao 3 

500 

10.3 

220 

310 

. 1 

80 

8.1 

100 

5000 

10.3 

2.200 

3100 

.3 

2510 

14.0 


instrument errors is assumed to be 10 arc sec, and the interpolation errors 
were calculated to an accuracy o£ about 0.1 m/s and 1.0 km. 

The total position error [root sum square of random and interpolation 
errors) is smaller at the center of the interval except when r = /2 and 
At = 5, and in no case is it much greater than 1 percent of the total position 
vector. The total velocity error is also satisfactory both at the ends and 
the center of the interval, but it is much smaller at the center, indicating 
that interpolation to the center of the interval is preferable. 

The errors in velocity which arise from the use of linear interpolation 
were also calculated at the center and end of the interval and found to exceed 
the corresponding error from Gibb's method by one or two orders of magnitude. 
Thus acceptable accuracy might be obtained near the center of the interval, 
but linear interpolation is not a generally satisfactory method for computing 
velocity. 

At the smaller values of r the interpolation intervals (2At) used in 
table 7 would be no greater than the time required for two sets of observa- 
tions; but for larger values of r the time interval becomes much greater. 

For example, when r = 10, At in table 7 is 50 minutes, and if At were 
reduced to 10 minutes Oy wou ld be about 50 m/s or about 2 percent of the 
total velocity. For r = /103 the random error for At = 10 minutes would 
exceed 0.5 km/s. The position errors at the larger ranges appear to be rather 
large, but the root sum square of random and interpolation errors is less than 
1 percent of the total range. 

Equivalent data were not obtained for orbits of higher eccentricity with 
the exception of one parabolic case. The interpolation errors in this case 
are somewhat higher than for the circular orbit, but a marked reduction in 
velocity errors is obtained by using large time intervals in this case also. 
For still higher eccentricities the orbit is essentially linear at long ranges 
and the velocity interpolation errors will be quite small, again indicating 
the desirability of large interpolation intervals. 
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The results just discussed indicate that all of the measurement types, 
except possibly type 4, can provide satisfactory preliminary trajectories, 
provided the time interval between the two position fixes is chosen proper] y. 
For r > v^. At should be 10 minutes or more and the extended schedule would 
be of little value. Also, the interpolation should be to a time near the 
center of the interval even though this choice reduced the velocity error at 
the expense of some increase in position error. If r > 10, the time interval 
required to reduce the random velocity error is so large that it would be best 
to try to reduce the position error with redundant data before attempting to 
determine velocity. The development of such a scheme, however, is beyond the 
scope of this study. 

For near-planet orbits, where r < 4, the interpolation error from Gibb's 
method is not excessive if the interpolation time is near the center of the 
interval. Linear interpolation would be the major source of error, and this 
error could probably be constrained within acceptable bounds. 


Near-Planet Orbits 

The data presented have been restricted to r ^ /Z and special 
consideration must be given to cases with r < /Z. In this case the subtense 
angle exceeds 90° and the planet disk fills the field of view considered in 
calculating data, thereby excluding all of the measurement types except type 3. 
While it is theoretically possible to measure angles larger than 90°, it is 
questionable whether it would be practical to measure subtense angles of much 
larger magnitude. 

Since type 3 measurements do not require a view of the complete disk they 
can be used at small values of r, but the angle B, separating the planes 
defined by r and the unit vectors to the two stars, must be less than the 
arc of the planet's disk lying in the primary field of view. The field of 
view of the trisextant could probably be increased substantially from the 3° 
of the instrument described in reference 7, and if sufficient accuracy could 
be obtained, type 3 measurements would be atractive for use in near-planet 
orbits. These measurements in effect combine the sextant and stadimeter 
measurements of reference 3, and it should be possible to develop a similar 
manual computation method for a limited range of eccentricities. The practi- 
cality of such a system depends, of course, upon the ease and accuracy of 
making the measurements for T</Z. 

CONCLUSIONS 


An error analysis has been carried out for determining preliminary 
trajectories using five different types of measurements. The measurement 
type (type 1) consisting of simultaneous measurements of the right ascension, 
declination and subtense angle of the planet has been used as a standard of 
comparison. The following conclusions can be drawn from the results of the 
study: 
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1. If the stars used in the measurements can be chosen properly, two of 
the remaining methods, the photographic method [type 2) and the trisextant 
(type 3) need a single measurement to provide a position fix having errors 
within a factor of two of the standard of comparison. 

2. Similar results can be obtained using a standard space sextant 
(type 5) to make seven angular measurements and interpolating to achieve a 
common epoch. 

3. The measurement type (type 4) which makes simultaneous measurements 
of one star- limb angle and the subtense angle at three different times and 
interpolates is less accurate than the others and offers little advantage. 

4. It is best to choose stars such that the projection of the radius 
vector f onto the plane of the two stars lies in or near the smaller angle 
between them. 

5. The proper choice of stars would be enhanced by increasing the 
maximum measurable angle from the present 55+° to 90°. 

6. For magnitudes of r greater than ^ planet radii and properly 
chosen stars, the random position error can be approximated roughly by the 
radial error from a single subtense angle measurement. The random velocity 
error is roughly the position error divided by half the time interval between 
the two position vectors used to determine the velocity. 

7. Bias errors, particularly those due to irradiance, can be very 
important and every effort should be made to reduce the unknown biases well 
below the standard deviation of the random errors. For example, an instrument 
with a standard deviation of 10 arc sec is generally regarded as practical and 
it would be desirable to reduce the unknown biases to 2.5 arc sec or less. 

8. For reasonably small bias errors at distances from /2 to 200 planet 
radii from the central body, satisfactory preliminary trajectories (position 
and velocity known within a few percent) can be found using any of the measure- 
ment types, provided suitable stars can be chosen. 

9. For all trajectories except near-planet orbits the hand-held space 
sextant can provide satisfactory measurements for determining a preliminary 
trajectory. Furthermore, for sufficiently large radial distances from the 
planet, linear interpolation is sufficient for determining position, but a 
higher order formula such as that from Gibb's method should be used for 
velocity. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. , 94035, July 9, 1971 
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APPENDIX A 


EQUATIONS FOR POSITION FIX 


It is pointed out in the text that some types of measurements require 
interpolation to refer the measured angles to a common epoch. The equations 
for the position fix are derived in this appendix under the assumption that 
the interpolation has already been carried out or, equivalently, that all of 
the angles used were measured simultaneously. 


DETERMINATION OF RANGE 


For all five types of measurements the radial distance, r, of the planet 
center from the spacecraft is ultimately determined from the semisubtense 
angle, S, as shown in sketch (d) . All the equations in the report involving 



Sketch (d) 


r have been normalized by defining the unit of length to be the radius of 
the observed planet, that is, by setting R equal to unity. Consequently, 

r = CSC S (Al) 

COMPUTATION OF THE UNIT VECTOR u 
Type 1 Measurements 

Type 1 measurements are direct measurements of RA, D, and S, and it is 
shown in the text that 


Ui = cos 

RA cos 

“1 


U 2 = sin 

RA cos 


(A2) 

U 3 = sin 

D 

) 
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Type 2, 3, 4, and 5 Measurements 


These types of measurements are direct measurements of the star- limb 
angles Yi and Y2> and^ with the exception of type 3, the direct measurement 


z 



The angle, <j>, between the stars is t 
shown by spherical trigonometry that 


of 2S. For type 3, B is measured 
instead of S, which must then be 
computed from Yij Yz> 
details of this computation will be 
discussed in the next section, but 
for the present it will be assumed 
that S is known. Figure 1 is 
reproduced for aid in deriving the 
equations for the components of u. 

First the star-center angles 
and 1(^2 are computed using the 
relationship 

i|ji = Yi + 

where = 1 for Yi measured to 

the limb nearest the star and kj^ = -1 
for measurements to the far limb, 
from tabulated data, and it can be 


cos ipi = cos RA cos D 
cos ip 2 = cos (4> - RA)cos D 


(A3) 


or 


cos ip 2 = cos RA cos D cos (j> + sin RA cos D sin <j) 


(A4) 


Substitution in (A3) and (A4) for cos RA cos D and sin RA cos D from 
equations (A2) gives 


ui 


U2 


cos ifji 

cos ll )2 - cos cos Ip 
sin (j) 


(AS) 

(A6) 


Since u is a unit vector 


U3 


± A 



(A7) 


Combining equations (AS), (A6) , and (A7) into a single vector equation gives 
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cos 

cos \p2 - cos ^2 COS <p 
sin c|) 

± /l - - U2^ 


CA8) 


Special Computation for Type 3 

In the case of type 3 measurements, yi j Y 2 > 3 -nd B are measured 
simultaneously, and S must be calculated. The law of cosines for spherical 
triangles gives 


cos (j) = cos cos ip2 + sin ipi sin 1^2 cos B 


(A 9 ) 


or 

cos <p = cos(yi + kiS)cos(Y2 + ^28) + sin(Yi + kiS)sin(Y2 + k2S)cos B 

It can be shown that since ki = ±1 

cos (Yi + k^S) = cos (k^Yi + S) 
sinCYi sin(kiYi + S) 

so that 


cos tf) = cos(k^Y]^ + S)cos(k2Y2 

+ k]^k2 sin(k^Yi + S)sin(k2Y2 S)cos B CAIO) 

Expanding the trigonometric functions of sums of angles in equation (AlO) 
with the appropriate formulas gives 


cos cj) = [cos kjY]^ cos k2Y2 kik2 sin k^Y^ sin k^Y2 cos B)cos^ S 

-[sin kjYi cos k2Y2 cos k;^Yi sin k2Y2) Cl ~ kik2 cos B)sin S cos S 
+ [sin kiYi sin k2Y2 cos k^Yi cos k2Y2 cos B)sin^ S [All) 


We note that 
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cos^ S = 


1 + cos 2S \ 


sin"^ S = 


sin S cos S = 


1 - cos 2S 
2 

sin 2S 


(A12) 


so that 


2 cos (j> = cps(kjY;^ - k2Y2) (1 ^1^2 ^os B) 


+ cos(k^Yi + ^2^2 ’’’ (1 - k^k2 cos B) 


cos(k^Yi + ^2"<'2 


2 cos (j) - cos(k]^Yi - ^2Y23 (1 + *^1^2 cos B) 
1 - kik2 cos B 


(A13) 


Once equation (A13) has been solved for S we can compute ipi and \p2 use 
equations (A5) and (A6) for u^ and U2- Equation (AS) can be used for U3 or 
the law of sines for spherical triangles can be used to show that 


sin D 


sin ijj^ sin t|^2 ‘S’ 

sin B 


so that 


U3 


= + 


sin ipj, sin sin 


sin B 


(A14) 


Equations (A5) , (A6) , and (A14) can be combined to give 


“ 


— 

cos 

'Pi 


cos 

'Pz 

- cos 4’l cos (j) 



sin (f> 

±sin 

'Pz 

sin \pi sin B 



sin (|) 


(A15) 


Equations (A8) and (A15) are equivalent and the one to be used for type 3 
measurements could be chosen for ease in computation. 
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APPENDIX B 


INTERPOLATION FORMULAS FOR POSITION AND 
VELOCITY FROM GIBB'S METHOD 


Following the approach given in reference 9 (p. 147), we expand the 
position vector, r, in a series in time as follows: 

r(t) = a + bt + ct^ + dt^ + . . . (Bl) 

We assume that r(ti) and r(t 3) have been calculated from measurements and 
wish to determine f(t2) by means of interpolation. For convenience assume 
that t2 = 0 and write 


r(ti) = a + bt^ + cti^ + dti^ 
r(t 2 ) = a 

r(t 3 ) = a + btg + ctg^ + dts^ 

I 

Differentiation of equation (Bl) twice and evaluation at 

gives 


(B2) 


^ 2 > ^3 


r(ti) = A(ti) = 2 c + 6 dti 
r(t 2 ) = A(t 2 ) = 2 c 
f (t 3 ) = Acts) = 2 c + 6 dt 3 > 


(B3) 


where A is the vector of gravitational force acting on the spacecraft, which 
can be calculated if the corresponding r (plus necessary ephemeris informa- 
tion in n-body cases) is known. Therefore A(t2) is_unknown, but the first 
and third of equations (B3) can be solved for c and d in terms of known 
quantities to give 


- ^ 1 rt3A(ti) - tiA(t3) 

“ 2 L t3 - ti 


^ ^ rA(t3) - A(ti) - 
6 L ts - ki 


(B4) 


and from equations (B2) 

[r(t 3 ) - r(ti)](t3 - ti) = b(t3 - ti) + c(t32 - ti^) + d(t33 - tjS) (B5) 
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Equations (B4) and (B5) can be solved for b, which is equal to the velocity 
at t 2 , giving 


b = v(t2) 


1 p 

t3"tl 


r(t3)-r(ti) 


A(ti) 


(2t32+2tit3-ti2) - 


A(t3) 


(t3^-2tit3-2ti2) 


(B6) 


We can solve for _a by substituting for b, c, and d in the equation for 
either r(ti) or r(t 3 ) giving 


a = r(t 2 ) 


1 

t3-fl 


t3i‘(ti)-tir(t3) 


A(tn 

6 


'tl't3 ('ti-2t3) 


tit3(2ti-t3) 


(B7) 


Note that equation (B6) can be solved with either t^ or t 3 set equal to t 2 , 
which has been taken to be zero; that is, given r(ti) and f(t 3 ) Gibb's method 
can be used to find the velocity at either end of the interval. In this case 
interpolation of f is unnecessary, and the solution of equation (B7) is 
trivial . 
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APPENDIX C 


FORMULAS FOR EVALUATING MEASUREMENT ERRORS 


This appendix is divided into four sections. The first section shows 
that the acceleration terms resulting from Gibb's interpolation method can be 
neglected in computing the errors in position and velocity. The next section 
gives derivations of the formulas used for evaluating position errors result- 
ing from measurement errors when only enough data are obtained for a single 
position fix. The third section presents derivations of equivalent formulas 
for the extended schedules. Finally equations are developed for evaluating 
the errors in determining velocity. 


SIMPLIFICATION OF ERROR FORMULAS FOR GIBB'S METHOD 


The purpose of this section is to show that the acceleration terms from 
Gibb's method can be neglected for the purposes of computing errors resulting 
from instrument inaccuracies. We differentiate equations (7) and (10) to get 

tsdr (ti)-tidr(t3) - tit3(ti-2t3) _ t2t3(2ti-t3) 

(Cl) 

and 




dv(t2) = 


dr(t 3 ) - dr(ti) 


dA(ti) 


(2t3^ + 2tit3 - ti2) 


dA(t3) 


(t3^ - 2tit3 - 2ti2) 


(C2) 


If conic orbits are assumed and 


r = ru 


then 


A = - 


yu 


and 


dA (tp) 


du(tj) _ 2u(tp)dr(ti) 

Lr2(t^) r3(t^) 


Substitution for dA(ti) and dr(ti) in (Cl) and (C2) gives 
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= Tr~rr |t 3 ^(tl)du(ti)(l + ^ ) + t 3 u(ti)dr(ti) (1 - Pi) 
-tir(t3)du(t3) (1 + - tiu(t3)dr(t3) (1 - P 2 ) 


dv(t2) = 


(■-») 

^ j^r(t3)da(t3)^l + 

('•^) 


) 


t3 

-r(ti)du(ti) 


r(t 3 )du(t 3 ) (1 + -Yj + dr(t 3 )u(t 3 ) (1 - P 3 ) 
dr(ti)uCti) (1 - P 4 ) 


) 


where 


3r3(ti) 


3r^(t3) 


P 

3r3 (t3) 


P 4 = 


3r^ (ti) 


(tl^ - 

2tit3) 


Ct3^ - 

2tlt3) 


(t3^ - 

2tit3 - 

2ti2) 

- 

2tit3 - 

2t32) 


We wish to know when the Pj^ are negligible with respect to unity. If we let 
ti = t 3 - 2At, thereby defining 2At as the interval between ti and t 3 . 


< 2At 

and 




|ti^ - 2 tit 3 | 

= [4At2 - 

t 3 ^| < 4At2 


|t3^ - 2 tit 3 | 

= j4t3At 

- t3^| < 4At2 

|t3^ - 

2tit3 - 2 ti 2 | 

= |- 8 At 2 

+ 12 t 3 At - 3t32| < 8 At 2 

|ti^ - 

2 tit 3 - 2 t 3 ^| 

= |4At2 - 

■ 3t32| < 8 At 2 

wish 

to know when 

8pAt^/3r3 

is negligible compared 


The values of p for the major bodies in the solar system, in units of 
(planet radii) ^/sec^ , are tabulated below in order of decreasing size. 
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Body 

Earth 

Mercury 

Venus 

Mars 

Moon 

Neptune 

Sun 

Jupiter 

Uranus 

Saturn 


p (planet radii) ^/sec^ 

1.54x10"6 

1.49x10~6 

1.37x10"6 

1. 18x10"6 
0.993x 10"6 
0.403x10-® 
0.394x10-® 
0.349x10-® 
0.280x10-® 
0.112x10-® 


This table shows that the maximum value of p is that of the Earth so that 


3r® 



If the acceleration term is to be less than 0.1, then, in seconds. 

At < 156 r®'2 

and for r = the interval (2At) between t^ and t 3 is about 8.4 minutes. 
This would require two type 3 measurements no more than 8.4 minutes apart 
while two independent sets of type 5 measurements would require the individual 
measurements to be about 1 minute apart. For r = 2, At is increased to 
7 minutes, requiring type 5 measurements at 2 minute intervals. If interpola- 
tion is carried out to the center of the interval, the maximum acceleration 
term is reduced to pAt^/r® and At is increased by 60 percent. Therefore, 
even for near planet orbits, it is reasonable to neglect the acceleration 
terms in most cases. Equations (Cl) and (C2) may be written 


dr(t2) 


t3dr(ti) - tidr(t3) 
^3 - ti 


(C3) 


dv(t2) 


dr(t 3 ) - dr(ti) 
^3 - 


(C4) 
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POSITION ERRORS FOR REGULAR SCHEDULES 


In the text the radial and normal components of the position error due to 
measurement bias are defined as br and bj^, respectively. To compute bj. and 
bj^ set the random measurement errors equal to zero in which case 


bj. = dr 


and (see eq. (13)) > 

bj^ = r/iui ^ + du2^ + du3^ j 

Since the random measurement errors produce no mean position errors, the 
standard deviations of the components of position error are given by 

= /e (dr2) - bj.2 I 

Oj^ = /r^E(dui^ + du2^ + du3^) - b^^^ | 


(C5) 


(C6) 


the equations for b^, bj^, a^. and are derived by substituting the measure- 
ment biases and random errors into the appropriate differentials as shown in 
the following paragraphs. 


Type 1 Measurements 
For this type of measurement 


r = CSC S 


(C7) 


and differentiation gives 

dr = ctn S esc S dS = -r/r^ - 1 dS (C8) 

From equations (14) in the text 


ds = 6s + bj + bg 


Therefore, 


E(dr2) = [OgZ + (bi + bs)2]r2(r2 - 1) 


and 

bj, = -(bj + bs)ri/r2 - 1 


so that 


a-r = 0(-r 




(C9) 


(CIO) 
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From equations (2) in the text 

ui = cos RA cos D 

U 2 = sin RA sin D 

U 3 = sin D 

Differentiating and substituting for dRA and dD gives 

dui = + bj^)sin RA cos D + (Sq + bp) cos RA sin D 

du 2 = ('Sra * tij^)cos RA cos D - ( 6 q + bp)sin RA sin D 
du 3 = ( 6 £) + bj^)cos D 

so that 

du^^ + du 2 ^ + du 3 ^ = ( 6 p + bp)^ + (<5j^ + cos^ D 

therefore 

bj^ = + b^ cos^ D (Cll) 

and 

On = r/ap2 + a2^ cos^ D (C12} 


Basic Equations for Types 2, 3, 4, and 5 

All the types of measurements except type 1 use the same basic position 
fix equations, namely (C7) and 


Ui = cos 


(C13) 

„ _ cos ipz 

02 

- cos \pl cos cj) 

(C14) 


cos cj) 


U 3 = ± A - 

2 2 
- U2^ 

(CIS) 


It is shown in appendix A that for type 3 we can also write 

sin \[ii sin ip 2 sin B 


but the derivative of this expression and that of (CIS) are identical so only 
(CIS) will be considered. 
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The expression for dr is given in equation (C8) and differentiating 
equations (C13) , (C14) , and (CIS) gives 


dui == -sin dipi 


du2 = 
du3 = 


sin ip 2 dip 2 


sin <p 
(u^duj + U2du2) 


dui ctn <f) 


From equations (CIS) and (C18) 


(C16) 

(C17) 

(C18) 


du^^ + du2^ + dug^ = — - — ^('1 _ j. 

u. 


[(1 - U2^)dui^ + (1 - ui^)du2^ + 2 uiU 2 dujL du2] 


Substituting for du^ and du2 from (C16) and (C17) gives 

dui^ + du2^ + du3^ = I d^il^i sin^ ’1^1 [ (1 ~ ^2^) + (1 - Ui^)ctn^ (|) 

U 3 * 

-2 uiU 2 ctn ([)] + dijj 2 ^ sin ip 2 ^ ^ 2j 

sin^ (p 

2dibi dip^ sin sin ^9 r r. . n) 

■ s'iiri^ ^ctn <p - U1U2]} 


If substitutions are made for u^ from (C13), U2 from (C14), U3 
(A14) and cos (p from (A9) , equation (C19) is reduced to 


(C19) 

from 


2 

^(dui^ + du2^ + du3^) = — (di(/i^ + dij;2^ - 2di/;i dip 2 cos B) (C20) 

sin^ B 


The appropriate expressions for dij;i, dijj2j and ds can be substituted into 
equations (C8) and (C20) in order to find bp, bn, Up and for the 
different measurement types. 


Type 2 Measurements 

In this case measure Yi, Y 2> ® simultaneously. The equation for 

dr is the same as for type 1 so that (C9) and (CIO) are used for bp and 

From the definition of the star-center angles 
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(C21) 


'I'l = Yi + kiS I 
'f'2 = Y2 ^^2^ j 

Differentiating and substituting the differentials from equations (12) gives 

dij^i = (6^ + bi) + ki (65 + bg) | 

d}p2 = (<52 + ^ 2 ) + ^2(<5s ^s) ) 

Substitution of dijji and d<j ;2 into (C20) gives 


(C22) 


bn = g [(bi + kibg)^ + (b2 + k2b3)2 - 2(bi + kibg) (b2 + k2bs)cos 

(C23) 

and 


n sin B ^ Y 


[2a 2 + 20s^(l - kik 2 cos B)] 


1/2 


(C24) 


Type 3 Measurements 

With this type of measurement yi> Y2> B are measured and S must 
be calculated from equation (6) which is 

COS(kiYl + k2Y2 + 2S) = ^ L- : ^^2Y jjJl,,l kj,k2_COS ^ 

(1 - kjk 2 cos B) 

If this expression is differentiated, it is found that 

rsin(k^Yi - k 2 Y 2 ) Cl + b^k 2 cos B) 


ds = 


1 

2 


(kidYi +k2dY2) + (kidYi - k2dY2) 


sin(kiYi+k2Y2'*‘2S) (l-kik 2 cos B) 


dBkj^k^ sin B[2 cos(k^Yi - k^Y 2 ) " 2 cos ([] 
sin(kiYi + k 2 Y 2 + 2S) (1 - kik 2 cos B)^ 

It can be shown from equations (C21) that 

*^lYi = ^i 4 ’i ■ ^ I 

k2Y2 = k2'j^2 - S j 

Substitution of equations (C26) into (AID) gives 


(C25) 


(C26) 


cos (j) = cos kn|;i cos k2^z kik 2 sin sin k2<J^2 B 


(C27) 
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Equations (C26) and (C27) can be used to reduce (C25) to 

. If.,. - k2i|;2Hl +kik2 cos B) "I 

ds = - 2 ((kl dYi +k2 dY2) - (ki dYi -k2 ^^3 g) J 

^ 2 sin i|;i sin sin B dB 1 

sin(knJ;i + k2'|^2) (1 ” kik 2 cos B) ) 

If we define 


and 


then 


P, = 


sin(kiij;i - k2^2^ "*■ kik 2 cos B) 

1 ~ sin(kii|)i + k2<j^2) (1 " kik 2 cos B) 


sin sin 4'2 ^ 

^2 sin(kit|;i + k 2 <J^ 2 ) (1 " kxk 2 cos B) 


ds = - j [ki dYi(l + p^) + k 2 dY 2 (l - P^) + 2 dBp^] 


(C28) 


and 


/dyj - k,k2 dYo\ 

= dYi + ki ds = ( 2 ) ~ Pi) “ kip2 dB 


( dYo - k,ko dYi \ 

— ^ -j (1 + Pi) - k 2 P 2 dB 

Substituting for dYi, dY2, and dB from equations ( 12 ) gives 
ds = - y [ki( 6 i + bi)(l + Pi) + k 2 (i 52 + b2) (1 - Pi) - 2 bj + 2 ( 6 g + bg)p2] 


d 4 ^i 


dip 2 = 


rc^i 

- kl’^2^2) ^ 

(b^ - k^k2b2) ■ 

(1 - Pi) - 

(6b + bg)kiP2 

- 

2 

- 

■(62 

- kik26i) + 

O^Z ~ kik2bj) 

(1 + Pi) - 

(6b + bg)k2P2 

> 

- 

2 

- 


(C29) 


When equations (C29) are used with (C8) and (C20) , it is found that 


bj. = -r/r^ - 1 y [kibi(l + pi) + k2b2(l - Pi) + 2bgP2 - 2bj] (C30) 
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(C31) 


= r/r2 _ 1 (I + + cTg^P2^] 


1/2 


2 ^ ^ - Pl^^cos B] 


+ 2bg(kibi - k 2 b 2 )PiP 2 (l " ^ik 2 cos B) 


+ 2bg^p2^(l - kik 2 cos B) 

si^B ^ ^ ^^1^2(1 - Pi^)cos B] 

+ 2ag2p2^(l - kik 2 cos B)| 


1/2 


(C32) 


(C33) 


Type 4 Measurements 

This method consists in measuring S and at times t^ and ts and 

interpolating to find cos \p, at t 2 while we measure S and Y 2 t 2 - 

Therefore d ^2 and ds are the same as for type 2 and br and are given by 
equations (C9) and (CIO) . 


Cos \pi (t 2 ) is given by equation (9) as 


cos 1/^1 (t2) 


1 

r(t3 - ti) 




6r^ 




2t3) 


cos 


'(^l (tl) 


pt 1 t o 

tl - - — ^ (t3 - 2ti) 
6r3^ 


cos ijJiCts) 


It was shown earlier in this appendix that the term involving y/r^ can be 
neglected in computing the errors. Thus, replacing r by 1/sin S 


cos i|;i(t 2 ) 


sin S(t 2 ) r t 3 cos jpi (tl) cos (tg) - 
t 3 " tl _ sin S(ti) ^ sin S(t 3 ) 


(C34) 


Differentiate (C34) and then, since relatively small time intervals are 
considered, it can be assumed that 'l^i(ti) = 'j^i(t 3 ) = i|ji and S(ti) = S(t 3 ) =S. 
This gives 
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(C 35 ) 


= -dS(t2)ctn ctn S + 


t3 dS (ti) - ti dS (t 
t3 - tl 


?) 


ctn ctn S 


^ tg diJ^iCt ij - tl difjiCtg) 
^3 - tl 


Since the biases are constants, equation (C 20 ) can be used* to write 
d’l^l (ti) = [6i(ti) + bi] + ki[65(ti) + bs] and from equation (12) 


dS(ti) = 6s (ti) + bs + bj 


Substitution of these expressions into equation (C 35 ) gives 

rt^6s(ti) - ti6s(to)"| 

d>l^l(t2) = [ ^ (ki + ctn i/^i ctn S) 


rt3'6i(ti) - ti6i(t3)~ 

-6s(t2)ctn tpi ctn S + + bi + kibg (C 36 ) 


From equation (C 22 ) 


d4>2(t2] = [62 (t2) + k26s(t2)] + (b2 + k2bs) (C 37 ) 


It can be seen that the bias terms in equation (C 36 ) are the same as those for 
di|^l in equations (C 22 ) so that b^ is given by (C 23 ) . Substituting equa- 
tions (C 36 ) and (C 37 ) into (C 20 ) gives 


(a 2 

r ^32 + t,2 n 

+ 

1 

Ct 3 - ti)2_ 

S 


ctn'^ ipl ctn“^ S + 1 


t 2 + t 2 
Ct3 - ti)-^ 


ctn S)^ + 2k2 ctn ipi ctn S cos B 


1/2 


(C 38 ) 


The quantity - ti)^) varies from a maximum of unity for ti or 

t3 equal to zero to a minimum of 0.5 for -ti = t3 = At. Since the latter 
would be more nearly true in practice, the results will not be greatly 
degraded by making this approximation, in which case 
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= 


n sin B 2 Y 


lo .2 . 




3(r"^ - l)ctn^ i|;i + 3 


+ 2 /r^ - 1 ctn (2k2 cos B + 


1 /2 


(C39) 


If tp 2 is interpolated instead of ipi , it replaces in equation (C39) 


Type 5 Measurements 
The schedule for this method is as follows: 


Time 


^2 

t3 

t4=0 


t6 

ty 

Angle 

measured 

S 

Yl 

S 

Y2 

s 

Tl 

s 


If the intervals between the tp 
interpolation 

dS(t2) 

dS(t4) 

dS(te) 


are assumed equal, then from linear 


dS(ti) + dS(t3) 
2 

dSCtg) + dSCts) 
2 

dSCts) + dS(ty) 
2 ~ 


Since r = esc S(t43, equations (12) and (C40) give 

r p^s Ct3) + <5s (t5) 


dr = -r /r^ 


+ bj + b3 


(C40) 


from which 

bj. = - (bj + bs)r /r^ - 1 

and 



(C41) 


(C42) 
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Equation (C20) can be used for this type of measurement if ti, t 2 , and 
t 3 are replaced with t 2 , and te, respectively. Then using equa- 

tions (C40) for the appropriate dS(tj^3 gives the equation 

di|^l(t 4 ) = ctn ctn S - dSCtg) - dS(ts3 + dS(t 7 )] 


ki 


[dS(ti) + dS(t 3 ) + dSCts) + dSCty)] 


^ j~ dy^ (t^) + dy^ (tg) j 

Substitution for the differentials from equations (12) gives 

d'l^i (t4) - - <Ss (t3) - 6s (ts) + 6s (ty) ] + -^ [6s (ti) + 6s (t3) 

r6i (to) + 6i (to) q 

+ 6s (t5) + 63 (ty)] + 2 J t>i + kibg 

Similarly, it can be shown that 




+ 



2 


+ b2 + k2bg 


The bias terms in these equations are the same as for types 2 and 4 so that 
bn is given by equation (C23) . Using di|;i (t 4 ) and di|; 2 (t 4 ) in equation (C20) 
gives 

On = [(i’^ - l)ctn2 4^1 + 2 k 2 (/r^ - 1 ctn \pi - ki)cos B + 3 

(C43) 



Position Errors for Extended Schedules 

It was pointed out earlier in the computation of instrument errors that 
the acceleration terms from Gibb's method may be neglected to give a formula 
of the form 


dr(t2) 


tg dr(t]^) - ti dr(tg) 
t3 - ti 


(C44) 


from which 


_ tg2|df(tg)2 - 2titg[dr(ti) ‘ dr (t g ) ] + ti2|dr(tg)|2 

(t3-U,^ 

If we assume that t 3 = -tj = At then 
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I 


(C45) 


+ |dr(t3)|2 

By definition 

r = ru 

dr = dr u + r du 

and u and du are orthogonal. Therefore we can write 

|dr(ti)|2 = dr^(ti) + r2(ti) |du(tj^) |2 (C46) 

and 


dr(ti) • dr(t 3 ) = drCtJdrCts) [u(ti) • u(t 3 )] + r (t 3 )dr(ti) [u (ti) • duCtj)] 

+ r(ti)dr(t 3 ) [u(t 3 l • du(ti]] + r(ti)r(t3) [du(ti) ■ du(t 3 )] 


(C47) 

If the change in true anomaly between t^ and is small, then the 
approximations may be made that 

u(ti) • u(t 3 ) = 1 

u(ti) • du(t3) = u(t3) • duCti] = 0 
r(ti) = rCta) = r(t 3 ) = r 

Then equation (C47) reduces to 

dr(ti) • dr(t 3 ) = dr(ti)dr(t 3 ) + r^ du(ti) • du(t 3 ) (C48) 


where the first term on the right-hand side comes from the radial components 
and the second from the normal components. 


Equations (C46) and (C48) can be used to separate equation (C45) into its 
radial and normal components to give 


dr2(t2) = 


dr[ti) + drCtg) 


-.2 


= -r.2 


r-=^(r^ - 1) 


dS(ti) + dSCtg) 


(C49) 


and 


r^ i du(t2) I ^ 


r^[|du(ti)|2 + 2du(ti) . du(t3) + |du(t3)P] 

4 


(C50) 
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The terms r2|du(ti)p and r2|du(t3)p can be found by substituting the 
proper time arguments into equation (C20) , and if equations (C16), (C17), and 
(Cl 8) are used with the proper time arguments, 

[du(ti) • du(t 3 )] = — ^ {dij;i (ti)di|;i (t 3 ) + di /^2 (t 3 ) 

sin*^ B 

- [di|Ji (ti)di|)2(t33 + di/^i (t3)di(;2(ti)]cos B} 

Making these substitutions into (C50) results in 

r2|du(t2)l^ = - — g — {[dt|^i(ti) + dt|;i(t3)]2 + [di|)2(ti) + di|«2(t3)]^ 

4 sin^ B 

-2[dtjii (tjdtj;2(ti) + dti^i (t3)dip2(t3) + di/^i (tj )di/;2 (ts) 

+ di^/i (t3)di|j2(ti)]cos B} (C51) 

Since biases are constants, it can be seen from equation (C49) that bj. is 
the same as for a single measurement of S and, hence, is given by (C9). 
Likewise, examination of equation (C36) shows that the biases in are not 

affected by interpolation so that bj^ is given by (C23) . We can evaluate 
Or and for types 4 and 5 with the extended schedule by setting the measure- 
ment biases to zero and taking expected values in equations (C49) and (C51). 


Type 4 Measurements 

The extended schedule for this type of measurement is 

Time t^ t 2 t 3 

Angles measured YjjS Y 2 »^ 

Angle computed ijii ij^z 

In this case interpolate to get (ta) > and use it with and S(t 2 ) to 

compute r(t 2 ). Then interpolate to get for ase in computing f(t 4 ). 

The assumptions made in deriving equations (C49) and (C51) require that 
t 4 = -t 2 = At/2 and ts = ti = 3At/2. 



Substituting the appropriate time arguments into equation (C49) gives 

o o o rdS(t,) + dSCt.)-|^ 

dr2(t3) = r2(r2 - 1) [ 2 


and since S(t 2 ) and S(t 4 ) are measured directly the random components are 
independent and 
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(C52) 


/r^ - 1 
Or = ^3^/— 2 

If we omit the biases and assume the values of tj^ given above, then 

r6s(ti) + 6s (t.)-] 

dipi (t2) = I 2 — J(ki + ctn ipi ctn S) 

T 6 , (t,) + 6 , (t.) 1 
- 6 sCt 2 )ctn ipi ctn S + J 

d'|^2('t2) = 62 (t2) + k26s(t2) 

dijji(t4) = 61 (t 4 ) + ki6s(t4) 

d 4 » 2 (t 4 ) = I 2 J (^^2 ■*■ otn \p2 ctn S) 

r6i (to) + 6 ^ (ts)1 

- 6 s(t 4 )ctn ip 2 ctn S + I I 

When these differentials are substituted into (C51) we can show that 

On = 2 sin B ><'1 + 3 ctn^ ^2 

+ 4 ctn ipj ctn ^2 00 s B) + /r^ - 1 (4ki + 6 k 2 cos B)ctn 

) 1/2 

+ /r^ - 1 ( 4 k 2 + 6 k 2 cos B)ctn ^2 + ^0] | (C53) 

Type 5 Measurements 

The extended schedule for this type of measurement is: 


Time 


^2 



^5 

te 

ty 

^8 

tg 

Angle measured 

S 

Yi 

S 

Y2 

S 

Yl 

s 

Y2 

S 


This schedule is the same as for type 4 except that in order to compute ipi 
associated with each measured Yi we must interpolate S using the measured 
values at_ adjacent times. Then following the procedure used for type four we 
compute r(t 4 ) and rCtg) and then interpolate to get r(t 5 ). 
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Substituting the appropriate time arguments into (C49) gives 


dr^t3) = - 1) . dSCt^) 

and taking account of linear interpolation 

~ | ^dS(t3) + dS(tg) ^ dS(tg) + dS(t<;^)J 

Assuming the biases to be zero and taking the expected value gives 


s^'/ 8 


5(r^ - 1] 


(C54) 


Similarly, by making the proper substitutions into equation (C51) we can show 
that 




2 sin B Y 


0 Q ^ 

^0 2 -j- r — 1 T9 r*‘f"Ti2 ill, 4- 9 r*1"-n2 


[ (r^ - 1) (2 ctn"^ ipl + 2 ctn^ ijj2 


+ ctn rpi ctn 1 P 2 cos B] + 5/r^ - 1 (kj ctn 1^2 ^2 ctn ijJi)cos B 

1/2 

- 15k]^k2 cos B + 20 } (CSS) 


VELOCITY ERRORS 


The errors in velocity are determined by differentiating equation (10). 

It was shown earlier in this appendix that, for the purpose of computing the 
position and velocity errors due to measurement errors, the acceleration terms 
in equation (10) may be omitted. The resulting equation is of the form 


dr(tg) - dr(ti) 

dv2 = T — — 

t3 - tl 


(C56) 


Since the bias errors are constants, it can be seen from (CS6) that there 
are no errors in velocity due to measurement biases. Therefore we will com- 
pute the random errors in velocity under the assumption that the biases are 
zero. Hence, 


E[|dv(t2) p] 


I |dr(t3)[^ - 2[dr(ti) • dr(t3)] +|dr(ti)|2 

\ Ct3 - ti)2 


(C57) 
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If r(ti) and r(t 3 ) are determined from two independent sets of measurements, 
which is the case except when the extended schedules are used, the random 
components of dr(ti) and dffts) are independent, and the expected value of 
the dot product is zero. 

From equation (13) 


|dr(ti)|2 = dr^Cti) + r^ (tj[) [dui^ (ti) + dua^ (ti) + du 32 (ti)] 
so that substitution from (14) , assuming bj. = bj^ = 0 into (C57) gives 

E[|dv(t2)|] = T — / 2(01-2 + cr 2) 

t3 - ti / 


If At is defined b_y 


then 


t3 “ t^ = 2At 


E[|dv(ta)|] = 



Similarly, if the acceleration terms are neglected, 

t3df(ti) - tidf(t3) 


dr(ta) = 


At 


(C58) 


and 


dr(ta)|2 = • d?(t3)] + ti2|d?(t3)|2 

4At2 


Therefore, for two independent sets of measurements and assuming bj. = b^^ = 0 


E[|dr (ta)!] 



(C60) 


From equations (CSS) and (C60) 


E[|dv(ta) I ] 


E[|dr(t^) I ] 

/ti2 + t32 


(C61) 


The quantity (ti2 + t 32 ) varies from a maximum of 4At2 for interpolation to 
one of the ends of the interval, where ti=ta=0 or t 3 =ta= 0 , toa 
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minimiun of 2At^ if t 3 = -ti = At. Therefore, for interpolation to the center 
of the interval 


E[|dv(t 2 )|] ^ j_ 
E[|df(t2)|] 


(C62) 


While for interpolation to the ends of the interval the ratio is reduced by 
'/2. Note that this change is caused by an increase in the standard deviation 
of the position error while that of the velocity remains the same. 


For the extended schedules assume that t 3 = -t^ = At so that 
equation (C59) reduces to (C45) , while (C57) is the same except that the sign 
of the dot product is reversed. Making this sign change in the derivations 
for the standard deviations of the position error gives for type 4 measurements 


a 


a 


I 

r 


n 


t^sr /r2 - 1 
At / 2 


2At sin B Y 


5a 2 - 1) (3 ctn^ + 3 ctn^ ijj2 


-4 ctn ipi ctn ip2 B) + /r^ - 1 (4ki - 6 k 2 cos B)ctn 



I- /r^ - 1 (4k2 - 6ki cos B) ctn il >2 + 10] > 


For type 5 measurements 


^ / Cr^ - 1) 
At / 8 


■'n 


2At sin BY 8 


5a„2 + — [(r^ - 1) (2 ctn^ \pi + 2 ctn^ i|)2 


-ctn ctn xp2 cos E) - s/r^ - 1 (k^ ctn ^^2 ctn i|; 2 )cos B 

+ 15k2k2 cos B + 20] 



where a^.' and Ojj' are the standard deviations of radial and normal errors in 
velocity. Numerical evaluation over the range of measured angles considered 
indicates that if t 3 = -ti = At then equation (C62) holds for the extended 
schedule within a factor of 2. 
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APPENDIX D 


EQUATIONS FOR ANALYSIS OF ERRORS DUE TO LINEAR INTERPOLATION 


The purpose of this appendix is to compare the magnitude of position 
errors arising from the linear interpolation of the range r or a star-limb 
angle y with the errors due to measurement inaccuracy. Since the interpola- 
tion errors are a function of the spacecraft trajectory, the state of the 
vehicle must be accounted for. This has been done by finding the minimum 
value of r for which the maximum interpolation errors are less than some 
fraction, M, of those due to measurement inaccuracy. An expression for the 
error was found in terms of the gravitational constant y, the true anomaly 
0, the orbital eccentricity e, the range r, and At which is half the 
interpolation interval. The largest value of y, in (planet radii) ^/sec^ , of 
the major bodies in the solar system was used, and 6 was chosen to maximize 
the error. This choice of 6 was made in order to simplify the calculations. 
The minimum value of r was then found for various values of eccentricity and 
At. 


The errors in range, which can be evaluated directly from a Taylor series 
expansion of r, are considered in the first section. Errors due to linear 
interpolation of y are dealt with in the second section by an indirect 
method because the Taylor series expansion of y contains a singularity. 


Interpolation of Range 

If r(t) is known at times t^ and ts, an appropriate value, r(t 2 ), can 
be obtained for r(t 2 ) under the assumption that r is a linear function of 
time. Then 


r(t2) 


(tg - t^jrCt^) + (t^ - tijrCtg) 
^3 - ti 


and if t 2 = 0 


r(t2) 


tgrCt^) - tj^rCtg) 
^3 - ti 


(Dl) 


Since r is not a linear function of time, we can estimate the error in 
r(t 2 ) by expanding r in a Taylor series about t 2 . This gives 


r(ti) 


r(t2) + tir(t2) 

t(t2) + t3f(t2) 


t- 2 ^3 

+ — r(t2) + r(t2) + 


— r(t2) + — r(t2) 


(D2) 
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Substituting (D2) into (Dl) gives 


r(t2) - r(t2) = r(t2) 


r(t2) 


+ ts) 
3: 


or, dropping the time argument. 


Ar 


-r 


■*^1^3 


- r 


^1^3 


(tl + t3) 


(D3) 


This is the error in interpolating with a first-order formula as compared with 
a third-order formula such as those derived by Gibb's method. It should give 
a good indication of the total error when (t 3 - t^) is small. Note that if 
t 2 lies at ti or t 3 or half way between t^ and t 3 (-t^ = t 3 = At), the 
third-order term disappears. Define Ar 2 and Ar 3 as the maximum values of 
the second- and third-order terms, respectively, in equation (D3) and deter- 
mine the minimum value of r for which |Ar 2 | ^ Ma^. In this case is the 

standard deviation of position error due to a single subtense angle measure- 
ment of 10 arc seconds standard deviation. Next find the minimum value of r 

for which |Ar 3 | is negligible compared to |Ar 2 | and hence compared to 

mS-X IT13.X 

the measurement error. If this latter condition is true, then the second- 
order term is a valid representation of the total interpolation error. If the 
trajectory is assumed conic 


r 


r 


h2 

m( 1 + e cos 0) 

ye sin 0 
h 


ye cos 0 


2y^e^ sin 0 cos 0 yeh sin 0 
^ “ hr3 r'^ 

or substituting for r in the second term from (D4) 

... _ y^e sin 0 


hr' 


(3e cos 0+1) 


(D4) 

(D5) 

CD6) 


(D7) 


where e is the orbital eccentricity, 0 is the true anomaly and h is the 
angular momentum. 

From equations (D3) and (D6) the second-order error is 

, 'tit 3 ye cos 0 

'''■2 = - -2 ?2 
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and tits is maximum for ti = ts = At so that 


I Aral ^ 


At^pe 

2r2 


(D8) 


The standard deviation of the rand om erro r in determining range from a single 
subtense angle measurement is ar/r^ - 1 where a is the standard deviation 
of the instrument error. For the ratio of the maximum of lAral to this 
quantity to be less than or equal to some number, M requires 



peAt^ 

2oM 


(D9) 


The minimum values of r from equation (D9) are plotted in figure 17 for 
M = 0.1 and 10 arc seconds, and are discussed in the section on results in the 
text. 


From equations (D3) and (D7) 


(3e cos 0 * 1) (DIO) 

0 


Since h = /prp (1 + e) and planet radii are the units of length 


h > /p(l + e) 


(Dll) 


and since 


3 

sin 0 (3e cos 6 + 1) = y e sin 20 + sin 6 
sin 0 (3e cos 0 + 1) < ^ ^ 


(D12) 


Therefore substituting from equations (Dll) and (D12) into (DIO) gives 


[Aral < 


tit3(ti + tg) 
6 


e(3e + 2) 
r^ 2/1 + e 


(D13) 


If the interval between ti and t 3 is 2At then |tit 3 (ti + t 3 ) | is maximum if 
ti = -At(l ± I// 3 ) so that 


+ t3) 

6 


< 


2At^ 

9/3 
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then 


Ar. 


< 


At^ H 

9/3 


3/2 


e(5e + 2) 
2/1 + e 


(D14) 


Dividing equation (D14) by (D8) gives the ratio of the maximum value of lAr 3 j 
to the maximum value of | At 2 | as 


^^3 max _ 2/^ 5e + 2 At 
I ^^2 I max " 9*/3 vTTT ^ 


5 

^ 4 
o 

QJ 3 
c 

a. 2 
c 

E I 



At, mm 


Figure 19.- Minimum radius for neglecting third-order term 
in analyzing error due to linear interpolation of radius. 


The third-order term can be consid- 
ered negligible if this ratio does 
not exceed 0.1, which requires 


^ 20/u 
^ 9/3 


3e -»• 2 


,/r 


+ e 


(D15) 


The minimum values of r for various 
values of e and At, using p for 
the Earth, are plotted in figure 19. 
All these values are considerably 
smaller than the corresponding ones 
for I Ar 2 I ^ Mor, thus indicating that 
the second-order error term is a 
valid approximation for the total 
interpolation error. 


Linear Interpolation of Star-Limb Angles 

In this section of the appendix sets of equations similar to (014) and 
(D15) are derived for the case where linear interpolation is used with a star- 
limb angle instead of r. For this purpose it is convenient to deal with the 
star-center angle ip at first. The error in position due to errors in ij;i 
and \p 2 is normal to the radius vector and can be shown from equation (C20) to 
be 


r/Suy 


+ du- 




/ 


For simplicity interpolation will be assumed for only one angle C'Pi or 4 ^ 2 ) 
a time and A\p will be the interpolation error. Then 
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(D16) 


cJd 


ui^ + 


du2^ + 


dug^ = 


r h\jj 
sin B 


Since actually y instead of ^ is interpolated, the "interpolated" value ^ 
of i|'(t 2 ) can be written 


ip = 


tgYCti) - t^yCtg) 


- ti 


+ kS(t2) 


(D17) 


Equation (D17) can be written 


ijj = - kS + kS (t2) 


where 


~ t3i|;(ti) - tnJ^Ctg) 

'P = rr-T-7E 


and 


S tgSCtp - t,S(t3) 


It can be seen that ip and S are the values of ij^(t 2 ) and S(t 2 ) obtained by 
linear interpolation and 


/\\P = ip - \pCt 2 J 


= [ip - ip(t2)] - k[S - S(t2)] 


is the difference of the errors in interpolating \p and S. Furthermore, 


AiJj] < [i|) - if)(t23 1 + |S - S(t2] 


(D18) 
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The second term can be evaluated directly by means of a Taylor series 
expansion, but the first term requires an indirect approach. For this purpose, 
define ip with different geometric equations from those used in the text as 
illustrated by sketch (e) . The X-Y plane is the orbital plane of the 


z 



Sketch (e) 


spacecraft with the origin at the center of the planet and X axis at peri- 
apsis. The star-spacecraft angle is the compliment of the star-center angle, 
ip, defined in the text, and the elevation angle of the star above the orbital 
plane is defined as (180° - ipo) • Note that ipQ is the minimum value attained 
by \p during a complete orbit. The angle 0 is the true anomaly of the 
spacecraft, while A9 is the change in true anomaly from the spacecraft's 
position to the projection of the star on the orbital plane. 

For spherical trigonometry 


cos (180° - >p) = cos (180° - i|/q)cos A6 


which can be reduced to 


cos Ip = cos cos A0 


(D19) 


we can write 




= 0 


di/j 

dA0 
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so that if the time interval is small enough for 0 to be assumed constant, 
then A0 can be used as the independent variable instead of time. The mini- 
mum value of r for which this approximation is valid is found by expanding 
A0 in a Taylor series to give 


At20_ 

A0 = A0q + ^ + • • • 


and requiring that the second-order term be less than 10 percent of the first 
order. 


If the trajectory is assumed conic 

h 4r (1 + e) 
0 - 2 - 2 


(D20) 


0 


2ye sin 0 

r3 


(D21) 


then setting rp equal to one planet radius and [sin 9| to unity gives 

I At^ (0/2) I ^ e 

ht0| 

The minimum value of r for this ratio to be less than one tenth is given by 


^ lOv^ eAt 
r ^ — 3^33;:^ 
/l + e 



Figure 20.- Minimum radius for assuming d constant. 


the resulting value of r^in is 
plotted in figure 20 for various 
values of r and e using a 
geocentric orbit. 

If r is assumed to be greater 
than the minimum values in figure 20, 
then A9 can be used as the indepen- 
dent variable for interpolating ip. 
Since di(;/dA0 is singular for \p = 0, 
it is impractical to use a Taylor 
series, and we make further simplifi- 
cations by considering the plot 
(obtained from eq. (D19)) of ip as a 
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function of A6, shown in figure 21. It can be seen from the two examples of 
linear interpolation in this figure that the maximum interpolation error 
occurs when the epoch t 2 is at the center of the time interval and 
A9(t2) = 0- Under these conditions i|;(t 2 ) = 



Figure 21 Variation of \p with Ad and i//q, with two examples of linear interpolation. 


'P = a-rtd ^3 = “^1 = Therefore from equation (D19) 

cos ifi = cos cos A6 


and since 


cos i('(t2) = cos \pQ 


then 

cos ijj - cos 4^(t2) = -cos 'I'qCI - cos A6) 
and from trigonometry 

+ '^'o^ ■ 'f'r 


cos \jj - cos = -2 sin . 2 


sin 


From equation (D19) ip ^ \po ^ that 

'ip + 


sin 




sin \p 


(D22) 


(D23) 
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Therefore, combining equations (D22) and (D23) gives 


2 sin 


(V") 


< 


ctn if'oCl “ cos A0) 


or 


^ ctn 4>o sin^ ^ (D24) 

For the fastest orbit considered (u = 1.54x10 I’p = 1 and e = 10) it 
can be shown that the change in true anomaly from perigee in 5 minutes is 
about 0.89 radian. Since A6/2 ^ 0.445 radian, it is reasonable to make the 
approximation that 



sin 


2 ~ 2 


and from figure 21 


then 


- ’^o^ A6 



^ - 4^0 


2 


these approximations reduce equation (D24) to 


I ^ - '('o I = I ^ (^2) i < ^ ctn iJJq 


and if r is less than the appropriate minimum value from figure 19 

1 - 4' (t 2 ) I ^ ctn ipo 

We can replace 0^ with h^/r^ and eliminate h by using equation (04) to 
get 

q 2 = y(l + e cos 0) 

r3 
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so that 




< 


u(l + e)At^ ctn tJ'q 


(D25) 


Now consider the second term in equation CD18) , [S - S(t 2 )l- 
definition 


S = 


tgSCtj) - tjSCtg) 

t3 ■ 


and we can expand S(ti) and*S(t 3 ) in a Taylor series about t 2 to give 


S - S(t2) 


titg 


.. tito(ti + tq) 

S(t2] + ^ -i ^ S(t2) + . 


It was shown in the section on the linear interpolation of r that 

titgCti + tg) ^ 2At^ 

^ ^ 9/3 

and 


^1^3 < At^ 
2 '' 2 


so that 


|S - S(t2)| lS(t2)| + |S(t2)| + . . . (D26) 

9/3 

We wish to be able to approximate |S - S(t 2 )| by the second-order term and 
therefore find the minimum value of r for which the maximum magnitude of the 
third-order tenii does not exceed 10 percent of the second-order one. By 
definition 


sin S = 1/r 


where the equation has been normalized by using planet radii as units of 
length. Differentiation gives 

S = -r sin S tan S 

S = -r sin S tan S - Sf sin S(1 + sec^ S) 
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or, substituting for S 


S = -r sin S tan S + r^ sin^ S tan S (1 + sec^ S) 


Differentiating and substituting for S gives 


S = -sin S tan S[r - 3rr sin S (1 + sec^ S) 
+r^ sin^ S (2 + sec^ S + 3 sec^^ S)] 


Substitution of the derivatives of r from equations (D5) , (D6) , and (D7) 
gives 


S = sin S tan S 


y^e sin 9 
hr^ 


(3e cos 0 + 1) 


3y^e"^ sin 6 cos 0 . „ 2 

— + sin S(1 + sec-^ S) 

hr^ 


2y^e^ sin^ 0 


sin^ S(2 + sec^ S + 3 sec^ S) 


and substitution for sin B and tan S using the definition of sin S gives 


S = 


y^e sin 0 


hr 


1 L 


3e cos 0 + 1 3e cos 0 2 

+ (1 + sec^ S) 


_ (2 + sec2 S + 4 sec'^ S) 

h2 

Substituting from (D4) for r inside the bracket and from (D12) for h gives 

I's’l ^ ^ ‘iPl r^/r2 - 1 [sin 0[1 + e cos 0(7 + 3 sec^ S] | 

(1 + e) 


+e2 cos2 0(6+3 sec^ S) - e^ sin*^ 0 (2 + sec^ S + 3 sec** S) | (D27) 


2 • 2 
c 1 n ^ 


Substitution for r, sin S, tan S, and h in the expression for S gives 
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vie 



(1 + e)r^yx^ 


1 


e sin^ 0(1+ sec^ S) - e cos^ 6 - cos 9 


or 


|S| < — j== |e(l + sec^ S) - e(2 + sec^ S)cos^ 0 - cos 0| (D28) 

(1 + e)r^/r^ - 1 

We can rewrite (D27) and (D28) as 


|s| ^ 

(1 + e)^^^ 

and 

|S| < - 1 (D30D 

where F and G are defined as the maximum magnitudes of the terms in the 
absolute value signs. In order for the maximum magnitude of the third-order 
term in equation (D26) to be less than 10 percent of that of the second-order 
term 


r > F 

9/3 A + e ^ 


At 


(D31) 


approximations were used to find 


The values of F and G depend on 0 and sec^ S as well as e, and successive 

their values for various values of e and At. 
The resulting values of are 

plotted in figure 22 for geocentric 
orbits as functions of At for 
various values of e. If r is 
larger than these values and those in 
figure 20, equations (D25) , (D26) , 
and (D30) combine to give 



lAi/ii 




At^ y(l + e) 


ctn 


At , min 


yeG 

(1 + e)r^i^r^ 


t] 


Figure 22.- Minimum radius for neglecting third-order term 
in analyzing error due to linear interpolation of 
subtense angle. 
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This expression can be used with (D 16 ) to show that the upper bound on the 
position error r due to the interpolation of y is given by 


^ r yAt^ (1 + e) . , , 

< — : 5- ^ ctn 1^0 

sin B 2r^ 


yAt^eG 


2(1 + e)r^t^r^ 




(D 32 ) 


This error is compared with for type 5 since it has been 

found (see the section on results in the text) that this quantity is usually 
smaller for type 5 than for type 4 . From equations (C 42 ) and (C 43 ) 


Co/ . * | . [(^^) 


ctn“^ 


k2/r^ - 1 ctn ijji + k^k2 


cos B + 


Note that it has been assumed that ay = Og - a and that if i(j 2 is inter- 
polated instead of 4^1, it would replace ijji in the equation and k^ and k2 
would be interchanged. From this equation it is seen that dropping the sub- 
script on ip, that 






r*^ + 2 


2 4 sin^ B 

This expression is of the form 


/r^ - 1 ctn If; ^ (r2 - Ijctn^ 
2 sin^ B 


4 sin*^ B 


(Oj.^ + ^ r^a2(Ci^ ctn^ if; + 2C2 ctn ip + €3^) 


and from equation (D 32 ) 


fj2 > r^(Di^ ctn^ if; + 2D1D3 ctn ip + 03^) 


Therefore if the maximum m agnitude o f interpolation error is to be less than 
M times the maximum of 


Di < MaCi 
D1D3 < 

D3 ^ MoC3 
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Substituting for the Di and Cj^ gives 


yAt2 (1 + el ^ 


(D33) 


’'^^^^!^<M2a2/r2 - 1 

rS/j.2 _ I 


yAt2eG 

2 (1 + e)r2/r2 - 1 


sin B ^ Ma 


/ 


r2 + 2 


4 sin2 B 


(D34) 


(D35) 


It can be shown that if the first condition is satisfied, the second will also 
be, and only the first and last conditions will be considered. 

The minimum values of r, using y for the Earth and 10 arc seconds for 
a, for the two conditions were computed for each case and the largest taken as 
^min- These values of Vj^in plotted in figure 18 for two different 

values of sin B. The value of rmin from equation (D35) goes to zero with 
e so that for the smaller values of e the value of r,^in from (D33) is the 
larger and the curves do not change with sin B. For the larger values of e 
the errors increase rapidly as sin B decreases, and for very small values of 
sin B the interpolation error would become very large compared to the measure- 
ment error. However, the measurement error is also large in this region and 
such measurements would be avoided. 

Note that the values of r,„pj^ in figure 18 are always greater than the 
corresponding values in figures 20 and 22, so that the second-order terms used 
in obtaining the curves of figure 18 are valid measures of the error. The 
curves in figure 18 are discussed further in the text. 
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